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1.  This  report  is  thi  Final  Report  on  the  contract.  It  covers 
research  done  on  the  radar  cross  section  of  turbulent  wakes  during 
the  twelve  month  period  1 Jul  75  to  30  Jun  7b.  The  objective  of  this 
work  is  the  generation  of  a computer  code  which  can  accurately 
predict  the  radar  cross  section  of  turbulent  reentr/  wakes  of  arbitrary 
diameter,  electron  density  and  strength  of  turbulence. 


2.  The  above  work  is  of  value  since  it  provides  a method  for 
computing  the  radar  cross  section  of  ballistic  missiles  and  decoys. 
Including  the  doppler  spectrum  of  the  return.  It  will  be  used  by  the 
USAF  to  evaluate  the  possibility  of  discriminating  between  the  reentry 
wakes  of  large  and  small  vehicles. 
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I.  INTRODUCTION 


The  subject  of  this  report  is  a computer  model  for  the  radar  backBcatter- 
lng  from  a re-entry  induced  ionized  wake.  The  report  gives  both  the  theoreti- 
cal derivations  and  the  FORTRAN  coding  for  the  model,  and  thus  it  serves  as 
both  a final  report  and  a user's  manual. 

The  model  is  based  on  a previously  existing  model1,  which  applied  the 
same  basic  technique  but  in  a simpler  way.  Improvements  in  both  theory  and 
method  of  application  have  been  made  during  this  study. 

The  basic  technique  of  this  model  1b  to  use  a distorted  wave  Born  approx- 
imation in  which  the  nonrandom  background  medium  is  represented  as  a radially 
nonuniform  cylindrical  plasma.  The  propagation  of  electromagnetic  (em)  wa”es 
in  this  background  medium  is  formulated  rigorously  and  solved  numerically  to  a 
high  degree  of  accuracy.  The  results  of  these  calculations  are  used  to  "dis- 
tort", or  weight,  the  Born  approximation  for  the  scattering  by  the  turbulent 
fluctuations.  Section  II  gives  a more  detailed  general  description  of  the 
model  and  Section  III  gives  the  full  details  in  a manner  which  relates  direct- 
ly to  the  computer  code.  Section  IV  discusses  the  results  of  Borne  calcula- 
tions intended  to  demonstrate  the  model's  agreement  with  published  data  and 
to  show  the  significance  of  certain  phenomena  included  in  the  model.  Section  V, 
which  concludes  the  report,  gives  recommendations  for  future  development. 


ll  P.E,  Bisbing,  Final  Report:  The  EM  Scattering  Model  for  Re-entry  Wakes. 

Volume  I.  The  Wake  Field  Penetration  Model  of  Radar  Scattering,  "General 
Electric  Company,  RESD,  Report  No.  75SDR2386-I,  December  1975. 
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II.  GENERAL  DESCRIPTION  OF  THE  MODEL 
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A.  Background 

This  section  is  intended  to  give  the  theoretical  basis  for  this  model. 

This  theory  fits  within  the  general  topic  of  the  interaction  of  waves  with 
random  media.  In  particular  it  relates  to  em  wave. scattering  by  turbulent 
plasms.  Inasmuch  as  this  report  is  not  intended  as  a survey,  we  refer  the 
reader  to  the  review  of  this  topic  by  Granatstein  and  Feinstein2.  Also,  the 
closely  related  topic  of  em  wave  propagation  through  random  media  has  been 
surveyed  recently  by  Fante3.  Thus  we  will  refer  only  to  works  which  have 
established  those  points  which  are  pertinent  to  the  approach  used  herein. 

The  most  basic  part  of  the  model  is  the  Born  approximation,  which  ap- 
pears to  have  been  originally  applied  by  PekeriB4  to  acoustic  scattering  in 
turbulent  media.  Its  first  application  to  radio  waves  was  by  Booker  and 
Gordon8  and  to  ionized  media,  by  Booker6.  Tatarski7  gave  it  a generalized 
form  in  terms  of  the  turbulence  spectrum  function.  The  name  "Born  approxima- 
tion" is  borrowed  from  quantum  theory  to  denote  single  scattering;  i,e.,  it 
is  assumed  that  waves  propagate  undistorted  in  the  medium.  Improvements  to 
this  approximation  have  accounted  for  various  aspects  of  the  medium  which 
tend  to  distort  the  wave  before  (and  after)  it  reaches  the  scatterer  in 
question.  Thus  the  name  "distorted  wave  Born  approximation". 

Refraction  in  the  mean  wake  plasma  is  a very  significant  source  of  dis- 
tortion for  em  waves  in  ionized  media  because  it  tends  to  bend  the  path  of 
propagation  away  from  the  interior  of  the  medium.  In  fact  when  the  line  of 
incidence  makes  a low  angle  to  the  wake  axis  this  effect  is  important  at 
electron  densities  which  are  sufficiently  low  that  the  undlstorted  Born  approx- 
imation would  apply  at  normal  incidence.  Thus  propagation  in  the  background 
medium  is  the  main  source  of  distortion  used  in  this  model.  Other  distorting 
effects  are  Included  in  ways  which  are  similar  to  previous  theories. 

6.  A Green's  Theorem 

In  mathematicul  terms  the  Born  approximation  is  the  first  term  in  a per- 
turbation-iteration series  called  the  Neumann  series.  A survey  by  Keller®  gives 
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2.  V.L.  Granatstein  and  D.L.  Feinstein,  "Multiple  Scattering  and  Transport  of 
Microwaves  in  Turbulent  Plasmas",  Advances  in  Electronics  and  Electron 
Physics,  Vol.  32,  pp.  311-381,  Academic  Press,  N.  Y.,  1973. 

3.  R.L.  Fante,  "Electromagnetic  Beam  Propagation  in  Turbulent  Media",  Proc. 
IEEE,  Vol.  63,  No.  12,  pp  1669-1692,  Dec.  1975. 

4.  C.L.  Pekerls,  "Note  on  the  Scattering  of  Radiation  in  an  Inhomogeneous 
Medium",  Physical  Review,  Vol.  71,  No.  4,  pp.  268-269,  15  February  1947. 

5.  H.G.  Booker  and  W.E.  Gordon,  "A  Theory  of  Radio  Scattering  in  the 
Troposphere",  Proc.  IRE,  Vol.  38,  pp  401-412,  April  1950. 

6.  H.G.  Booker,  "Radio  Scattering  in  the  Lower  Ionosphere",  J.  Geophysical 
Research,  Vol.  64,  No.  12,  2164-2177,  December  1959. 

7.  V.I.  Tatarski,  Wave  Propagation  in  a Turbulent  Medium,  translated  by 
R.A.  Silverman,  McGraw-Hi]l,  N.Y.,  1961  (especially  Chapter  4). 

8.  J.B.  Keller,  "A  Survey  of  the  Theory  of  Wave  Propagation  in  Continuous 
Random  Media",  Proceedings  of  the  Symposium  on  Turbulence  of  Fluids  and 
Plasmas,  New  York-1968,  Polytechnic  Press,  Brooklyn,  1969  (pp.  131-142). 
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the  derivation  of  this  series  in  very  Bimple  terms.  Letting  M and  V represent 
the  nonrandom  and  purely  random  differential  operators,  respectively,  p the 
unknown  wave  function  and  g the  known  source  function,  the  equation  to  be 
solved  is 


(M  + V)p  - g . (1) 

Operate  on  both  sides  by  the  inverse  of  M. 

(1  + M“lV)u  - M~lg  (2) 

This  equation  can  be  written  as  follows: 

M - -M"lVy  + M“lg  (3) 

Returning  to  (2),  multiply  both  sides  by  the  inverse  of  the  binomial,  using  the 
binomial  theorem  on  this  inverse,  giving 

U - E (*-M“lV)nM“1g  (4) 

n*o 


which  is  the  Neumann  series . This  result  can  be  written  more  explicitly  as 
follows : 

U - E pn  (5) 

n-o 

y0  ■ M"lg  (6) 

yn  - -M-lVyn_1,  n>o  (7) 

where  pn  represents  the  order  contribution  to  y. 


Now  let  us  specify  the  above  system  of  equations  for  em  waves.  Assume 
all  quantities  have  the  time  dependence,  exp(-iut),  where  i is  the  imaginary 
number,  w is  the  radian  frequency  and  t is  the  time.  Then  the  vector  electric 
field  E obeys  the  differential  equation 


VxVxe  - k2(K+K')E  + ikJ 


(8) 


where  k is  the  free-space  wave  number  (2it  divided  by  the  wavelength),  K and  K' 
respectively  are  the  mean  and  random  components  of  the  complex  dielectric 
constant  of  tha  plasma  and  J is  the  source  (transmitter)  current  density  dis- 
tribution. This  equation  is  the  analog  of  (1),  where  y is  E,  M is  k2K-VxVx, 

V is  kzK'  and  g is  -ikJ.  The  inverse  of  M is  customarily  expressed  in  terms 
of  a dyadic  Green's  function  G,  which  satisfies  Maxwell's  equations  for  the 
electric  field  of  a unit  point  vector  source  in  the  mean  medium  only. 


VxVxg  - k2KG  + 4irlS(r-r' ) 


(9) 


In  this  equation  I is  the  identity  dyadic,  6 is  the  Dirac  delta  function  and 
r and  r'  respectively  are  the  radius  vectors  of  the  field  and  source  pointB. 
Apply  the  divergence  theorem  of  Gauss  to  the  combination  [Gx(7xe)-Ex(  Vxg)]. 
Use  of  certain  vector  identities  gives 


/ [G* (Vx7*E)  - E*(VxVxG)]dV  - / (VxE) • (GXdS) 

“/(VxG)*(EXdS)  (10) 

where  dV  Is  the  element  of  volume  enclosed  within  the  surface  having  outward 
normal  element  dS.  Let  this  surface  be  an  infinitely  large  sphere  containing 
both  the  source  and  the  scatterer  relatively  near  its  origin.  Then  both  E 
and  G are  outgoing  radiation  on  the  surface  and  the  curl  is  proportioned  to 
the  vector  product  with  the  normal,  so  the  right  hand  side  of  (10)  vanishes. 


Substitution  of  (8)  and  <9)  in  (10)  gives 

4ttE  - k2/K'E*GdV  + ik/J*GdV  (11) 

which  is  the  analog  of  (3),  so  M~l  is  an  integral  over  a Bcalar  product  with  G. 
The  analog  of  the  Neumann  series  is  the  following: 

E - Z En  (12) 

n«o 

4ttE0  - ik/J*GdV  (13) 

4^  - k2/K'Eft^1*GdV,  n>o  (14) 


In  summary,  (11)  is  a Green's  theorem,  in  which  G represents  propagation 
in  the  background  medium  (the  mean  or  nonrandom  component  of  the  total  medium). 
It  leads  to  the  perturbation  - iteration  solution  represented  by  (12)  through 
(14),  provided  this  series  converges. 

C.  Representation  of  the  Background  Medium 

Almost  all  previous  theories  ignore  the  background  medium.)  i.e.,  K is 
assumed  to  be  unity  everywhere.  Then  G is  the  free-space  Green's  function  and 
the  usual  Born  series  comes  about.  Accounting  for  the  background  medium 
obviously  is  not  easy,  but  in  the  present  problem  we  can  at  least  represent  it 
such  that  G can  be  calculated.  Thus  let  thu  mean  wake  be  an  infinitely  long 
cylinder  of  plasma  in  which  the  electron  density  and  collision  frequency  depend 
only  on  the  radial  distance  from  the  axis.  (The  effect  of  axial  gradients 
will  be  evaluated  later.) 

Consider  (13),  the  equation  for  the  aero-order  approximation.  If  the 
background  medium  were  free  space,  E0  would  be  the  field  radiated  by  the  trans- 
mitter because  G is  the  free-space  Green's  function.  The  effect  on  G of  adding 
the  inhomogeneous  background  medium  is  to  add  a scattered  wave,  which  also 
couples  to  J in  the  integral.  But  if  the  transmitter  is  in  the  far  field,  the 
coupling  of  the  scattered  part  of  G to  J can  be  neglected.  Thus  the  signifi- 
cance of  the  right  hand  side  of  (13)  is  that  E0  is  the  vector  field  induced  in 
the  background  medium  by  the  transmitter.  Note,  however,  that  it  is  the  total 
Induced  field,  incident  plus  scattered,  even  though  the  latter  component  is 
assumed  not  to  couple  to  the  transmitter. 

The  dyadic  nauure  of  G is  more  essential  in  (14),  since,  when  the  pertur- 
bation-iteration scheme  is  carried  out,  En  must  be  evaluated  locally,  at  all 
points  in  the  medium.  The  dyadic  qualities  of  G can  be  obviated,  however,  when 
the  n^  order  far  scattered  field  is  evaluated.  Simply  take  the  dot  product  of 
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both  sides  of  (14)  with  the  unit  vector  representing  the  receiver  polarization. 

This  means  that  En  on  the  left  represents  this  approximation  to  the  field  which 
*'  couples  to  the  receiver  and  that  G in  the  integral  is  the  total  (incident  plus 

scattered)  vector  electric  field  induced  in  the  background  medium  by  the  re- 
ceiver if  it  were  operated  as  a transmitter-  This  last  result  comes  about  be- 
cause of  reciprocity;  i.e.,  the  coupling  between  the  receiver  and  a small  test 
antenna  in  the  background  medium  is  independent  of  which  antenna  is  used  to 
transmit.  Therefore  the  effect  of  the  background  medium  is  quite  simple  in  the 
;•  first-order  approximation  of  the  far  field  since  one  needs  only  to  calculate 

( « the  vector  fields  induced  by  incident  plane  waves. 

!■;  _ In  order  to  lead  up  to  the  radar  cross  section  calculation  let  us  put  (14) 

!,  in  the  form  of  the  first-order  approximation  to  the  fields  scattered  by  an 

• incident,  unit  amplitude,  plane  wave.  Then  E0  in  the  Integra  is  the  total 

vector  elect  ic  field  induced  in  the  mean  medium  by  this  im  jnt  plane  y ave. 

The  fact  that  G is  the  fields  induced  by  radiation  from  the  i ceiver  implies 
f that 

4irrE  - k2eikr  /K'E0*G0dV  (15) 


In  this  equation  r is  the  distance  from  the  receiver  to  the  origin, of  the 
scatterer  and  U0  is  the  total  electric  field  induced  in  the  mean  toedium  by  a 
unit  amplitude  plane  wave  propagating  away  from  the  receiver  and  having  its 
phase  reference  at  the  origin  of  the  scatterer.  The  incident  waves  for  E0  and 
G0  have  the  polarization  of  the  transmitter  and  the  receiver,  respectively. 

We  now  give  the  problem  of  calculating  Eo  and  G0  precise  definition.  The 
incident  wave  propagates  at  an  angle  u to  the  wake  axis,  as  shown  in  Figure 
1.  The  axis  is  also  the  z-axis  of  a right  handed  cylindrical  coordinate  syBtem 
(p,  c|),  z),  where  the  direction  <j>  » o lies  in  the  plane  of  incidence  and  points  \ 

in  the  forward  direction.  The  complex  amplitude  of  the  incident  wave  has  the 
following  form:  (The  time  variation  will  be  suppressed  in  all  quantities.) 

I2(  =»  exp  [ik(p  sin  ix  cos  ^ + z cos  <»)  1 (16) 

1 i 

i 

Arbitrary  polarization  of  the  receiver  or  transmitter  is  represented  by  super- 
position of  two  principal  plane  polarizations.  The  electric  field  of  the  ; 

incident  wave  for  perpendicular  (often  called  transverse  below)  polarization 
is  taken,  by  arbitrary  choice,  to  point  into  the  plane  of  Figure  1 when  the  • 

wave  is  at  the  origin.  Thus  , ' 

Eit  “ iSi(P  sin  <|)  + $ cos  4))  (17)  ' 

where  12  ^ denotes  the  vector  incident  field  for  transverse  polarization  and 
denotes  ta  unit  vector.  Let  the  direction  of  parallel  polarization  be  such  that  ■ 

at  normn]  incidence  the  electric  vector  at  the  origin  would  be  aligned  with  the 

positive  z-axis. 

Eip  " Eilz  sin  a - (p  cos  4>  _ $ sin  4>)cos  a]  (18)  : 

With  the  assumption  that  the  electron  density  and  collision  frequency  depend 
only  on  p,  the  problem  of  the  background  medium  is  completely  specified.  Its 
solution  is  worked  out  in  detail  below. 
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D.  Representation  of  Turbulence  Effects 


Once  we  have  calculated  the  complex  vectors  K0  and  0o  at  all  points  In  the 
wake,  the  first -order  approximation  to  the  turbulence  scattering  la  given  by 
(15).  However,  a number  of  difficulties  are.  encountered  at  this  point. 

Perhaps  the  greatest  difficulty  with  (15)  is  the  statistics  since  K1  is  a 
random  variable.  The  most  direct  approach,  if  the  statistical  distribution  of 
K’  were  known,  would  be  a Monte  Carlo  analysis  in  which  the  fixed  distribution 
of  E0*G0  would  be  multiplied  by  each  realisation  of  K* . PerhapB  such  an 
approach  would  be  feasible,  but  it  has  not  been  attempted  in  this  Btudy.  In- 
stead we  have  relied  on  the  correspondence  between  (15)  and  the  Born  approxima- 
tion, which  occurs  when  the  background  medium  vanishes.  Thus  we  assume  that 
the  statistics  can  be  handled  in  a manner  similar  to  the  Born  approximation. 
Except  for  a spreading  owing  to  the  finite  volume,  the  Born  solution  is  propor- 
tional to  the  three-dimensional  Fourier  transform  of  the  turbulence  autocorrela- 
tion function  at  a wavenumber  which  belongs  to  the  vector  resultant  of  the 
product  E0*G0,  each  member  of  which  is  a unit  amplitude  plane  wave.  This  wave- 
number  is  ih  fact  equal  to  the  gradient  of  the  argument*  of  the  dot  product  of 
the  two  plane  waves.  Therefore  we  approximate  the  statistical  effect  of  turbu- 
lence by  taking  the  spectrum  function  at  the  wavenumber  which  is  equal  to  the 
gradient  of  the  argument  of  E0*G0. 

The  Born  result  is  proportional  also  to  the  volume  integral  of  the  mean 
square  fluctuation  of  the  dielectric  constant.  It  iB  obvious,  however,  that  if 
either  the  incident  wave  or  the  scattered  wave  were  attenuated  in  the  medium, 
this  volume  integral  would  have  to  be  appropriately  weighted.  Suppose,  for  a 
moment,  that  E0  and  G0  have  the  same  phaBe  as  the  incident  wave  but  different 
amplitudes  at  various  points  in  the  medium.  Than  the  expectation  of  the  radar 
cross  section,  which  is  proportional  to  tho  mean  square  of  (15),  would  be 
equal  to  the  Born  result  with  the  volume  integral  weighted  by  |E0*G0|2.  There- 
fore we  include  this  weighting  in  our  model.  Of  course  the  phaBe  of  Eo’^o 
varies  throughout  the  medium,  so  the  spectrum  function  at  the  appropriate  wave- 
number  weights  the  volume  integral  also,  rather  than  factoring  out  as  in  the 
Born  approximation  for  homogeneous  turbulence.  The  result  to  this  point  may 
be  written  in  terms  of  radar  cross  section  o as  Follows: 


C N 1 2 | E0  * G0 1 2 £ 
Att  r 2 1 

J 1 + (v/u>)' 


S (ke)dV 


S(k)S/D(r)eik*rdV  (20) 

where  re  is  the  classical  electron  radius  (2.81777  X 10“  ^m),  N'  is  the  electron 
density  fluctuation  (overbar  denotes  mean),  and  V is  the  electron  collison  fre- 
quency. In  the  definition  (20)  for  the  spectrum  function  S(k),  k is  a vector 
wavenumber,  r is  the  radius  vector  and  D(r)  is  the  autocorrelation  function  of 
the  electron  density  fluctuations.  The  effective  vector  wavenumber  ke  is  the 
gradient  of  the  argument  (phase)  of  Eq'Gq. 

Application  of  (19)  is  not  done  without  difficulty  in  certain  instances 

* The  term  ''argument"  denotes  thu  mathematical • term  for  a component  of  a 
complex  variable;  i.e.,  it  is  the  angle  in  radians  between  the  real 
axis  and  the  complex  vector. 
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owing  to  the  plasma  resonance  effect,  in  which  the  electric  field  tends  to 
infinity  at  the  surface  where  K tends  to  zero.  Thus  the  term  |eo*G0|z  causes 
the  integral  to  diverge  when  the  mean  wake  is  overdense  (electron  density 
greater  than  it/ (A ' rt.) , where  A in  the  free-apaoe  wavelength)  and  stiff  ie.iently 
eo!  1 is  Lon  Less.  In  the  prior  study  (see  Kef.  1)  we  did  not  ne  count  for  kt,, 

using  instead  the  argument  that  at  low  aspect  angle  the  wavenumber  was  not 
greatly  perturbed  from,  that  for  free  space  so  the  Born  spectrum  function 
factors  out.  Unfortunately  many  of  the  numerical  results  were  several  orders 
of  magnitude  greater  than  the  Born  approximation.  This  was  especially  true  in 
cases  for  which  the  diffraction  effect  allows  the  electric  field  to  penetrate 
where  refraction  would  otherwise  stop  it.  One  avenue  of  investigation  was 
followed  in  that  study  in  great  detail.  This  was  based  on  the  fact  that  the 
magnetic  field  and  the  product  of  the  dielectric  constant  and  the  electric  field 
are  not  resonant  even  when  the  electric  field  is.  Going  back  to  the  derivation 
of  a Green's  theorem  (see  Section  IIA),  the  analogs  of  (8)  and  (9)  can  be  writ- 
ten for  magnetic  fields  and  (10)  can  be  written  to  include  electric  and/or 
magnetic  fields.  Also,  (10)  can  include  K explicitly  by  applying  the  divergence 
theorem  to  Kn  times  the  type  of  combination  originally  used  above  to  derive 

(10) ,  where  n is  an  arbitrary  numerical  power  to  which  K is  raised.  One  can 
derive  an  unlimited  number  of  Green's  theorems  which  have  various  combinations 
of  K and  either  electric  or  magnetic  field  quantities  in  place  of  Ii  and  G in 

(11) .  Although  all  but  the  purely  electric  field  theorem  (11)  have  additional 
terms  involving  K and  sometimes  VK',  the  term  in  K'  is  always  similar  to  the 
first  term  on  the  right  hand  side  of  (11).  In  addition  to  the  result  of  (11) 
we  studied  the  results  of  theorems  using  H*F,  KE*F,  and  KE*KG,  where  H and  F 
are  the  magnetic  field  and  its  Green's  function,  respectively.  The  above-men- 
tioned additional  terms  were  neglected  and  radar  cross  section  was  calculated 
parametrically  for  each  of  the  four  theorems.  The  conclusion  was  reached  that 
negligible  benefit  could  be  obtained  from  these  alternatives.  In  all  cases 
where  the  resonance  effect  was  absent  for  any  reason,  all  four  theorems  gave 
essentially  the  same  results.*  In  other  cases  even  most  of  the  nonresonant 
theorems  tended  to  diverge  because  of  other  effects,  and  none  of  the  theorems 
agreed  on  the  answer. 

The  resonance  effect  has  been  used  by  Ruquist3  to  derive  a surface  scatter- 
ing effect.  But  we  find  a few  difficulties  in  that  approach  also.  In  the  dif- 
fraction regime  the  resonance  is  diffuse  in  that  the  electric  field  is  high  at 
not  only  the  critical  surface  but  also  over  a large  fraction  of  the  volume.  In 
the  geometric  optics  regime  the  critical  surface  is  inside  the  surface  where 
refraction  causes  ray  turning,  so  one  would  not  normally  associate  surface 
scattering  with  the  critical  surface  except  at  normal  incidence. 

The  approach  which  we  have  adopted  to  avoid  the  resonance  effect  is  to 
neglect  the  contribution  of  the  critical  surface  to  the  integral  in  (19)  at 
all  points  where  IEq^Gq]2  is  greater  than  unity  on  that  surface,  The  ration- 
ale for  this  is  based  on  the  fact  that  the  electric  field  has  a phase  reversal 
at  the  resonant  point;  i.e.,  the  phase  varies  infinitely  rapidly  through  the 


i 


* tTig  exception  to  this  conclusion  is  In  the  cross  polarized  backscatterlng , 
which  was  generally  considerably  smaller  when  the  magnetic  field  was  used. 

9.  R.l).  Ruquist,  "Resonant  Scattering  from  Weak  Fluctuations  in  Plane 

Stratified  Dielectric  Media",  Proc.  IEEE,  Vol.  63,  No.  1,  pp.  205-206, 
Jan.  1975. 
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critical  surface.  As  a result  the  wavenumber  goes  to  infinity,  where  the 
spectrum  function  vanishes  exponentially.  This  argument  is  admittedly  u bit 
weak  since  the  phase  of  E0*G0  does  not  have  a discontinuity  at  the  critical 
surface,  because  the  product  of  two  functions  which  chcnge  sign  at  a point  does 
not  itself  change  sign.  On  the  other  hand  it  is  found  In  numerical  studies  that 
the  effective  wavenumber,  defined  as  the  gradient  of  the  phase  of  E0‘G0,  does 
increase  rapidly  in  the  neighborhood  of  the  critical  surface.  An  example  of 
this  result  is  given  below. * It  is  interesting  to  note  that  in  some  theoretical 
studies  appearing  in  the  literature  the  effective  wavenumber  is  taken  as  the 
free-space  wavenumber  times  the  square  root  of  the  dielectric  constant  of  the 
background  medium.  In  that  approximation  one  would  have  zero  wavenumber  at  the 
critical  surface,  which  would  not  help  cancel  the  resonance  effect.  Of  course 
the  literature  which  uses  that  approximation  also  uses  geometric  optics,  which 
is  bothered  not  by  resonance  but  by  caustics. 

It  goes  without  saying  that  the  second  and  higher  order  approximations  for 
the  turbulence  scattering  are  very  complicated.  Thus  at  this  time  we  have  not 
developed  the  technique  for  rigorously  applying  the  model  to  the  multiple 
scattering  problem.  However,  in  order  to  include  the  main  effects  of  higher 
order  terms,  we  have  opted  for  an  approximate  representation  within  the  frame- 
work of  the  first-order  theory.  Two  phenomena,  attenuation  owing  to  random 
scattering  and  depolarization  from  multiple  random  scattering,  are  included  in 
this  representation.  The  equations  for  both  of  these  effects  are  taken  directly 
from  the  literature,  especially  in  the  case  of  depolarization.  A small  excep- 
tion is  the  fact  that  scattering  attenuation  is  not  Integrated  along  ray  paths, 
since  rays  are  not  used  in  this  model.  The  local  scattering  attenuation  is  used 
instead  to  define  a plasma  having  equivalent  attenuation  owing  to  absorption. 
This  plasma  is  then  linearly  added  to  the  original  mean  plasma  distribution  in 
order  to  simulate  the  effect  of  scattering  attenuation  on  the  coherent  wave. 


The  doppler  frequency  spectrum  of  the  scattering  1b  implicitly  represented 
by  (15)  when  the  theories  of  Lane10  and  Schapker11  are  invoked.  These  authors 
showed  that  the  spectrum  of  the  electron  velocity  fluctuations  in  the  wake  Is 
Gaussian,  which  Implies  that  the  first  and  Becond  moments  are  sufficient  to  des- 
cribe it.  Also  in  the  Born  approximation  the  doppler  spectrum  is  a linear  com- 
bination of  the  electron  velocity  spectra  in  the  wake,  so  each  of  the  doppler 
moments  is  the  same  linear  combination  of  the  corresponding  moments  of  the 
velocity  distribution.  ThuB  the  generalized  form  of  (19)  is 
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(21) 


where  Sn  is  the  moment  of  the  doppler  spectrum  and  un  is  the  moment  of 
the  electron  velocity  spectrum  at  a given  point  in  the  wake.  Of  course  the  radar 
cross  section  0,  as  given  by  (19),  is  the  zeroth  moment  (uq=1). 


*See  Section  111  C,  top  of  page  42. 

10.  F.  Lane,  "Frequency  Effects  in  the  Radar  Return  from  Turbulent,  Weakly 
Ionized  Missile  Wakes",  AIAA  J.,  Vol,  5,  No.  12,  pp  2193-2200,  December 
1967. 

11.  R.L.  Schapker,  "Frequency  Effects  of  Electromagnetic  Scattering  from 
Underdense  Turbulent  Plasmas",  AIAA  J.,  Vol.  8,  No.  9,  pp.  1582-1590, 
September  1970. 


E.  Representation  of  Axial  Variations 

Up  to  this  point  the  wake  has  been  treated  as  an  infinite  cylinder,  so 
provision  must  be  made  for  axial  variations  in  its  properties.  This  is  done 
simply  by  breaking  the  wake  up  into  a series  of  axial  stations,  each  of  which 
is  assumed  to  fit  the  uniform-  cyfinder  representation.  When  the  scattering 
from  any  given  station  is  calculated,  we  carry  out  the  volume  integration  in 
(21)  over  only  the  radial  and  angular  coordinates.  The  result  is  the  quantity 
dSn/dz,  which  represents,  for  example  if  n ■ o,  the  radar  cross  section  per 
unit  length  of  wake. 

An  integration  must  be  performed  over  the  axial  coordinate  z to  get  the 
total  scattering.  In  doing  this  we  assume  that  the  radar  signal  is  pulsed  and 
that  the  entire  wake  lies  within  the  main  beam  of  the  radar  antenna.  Then  the 
scattering  is  a function  of  the  apparent  wake  station  (time  delay) , which  is 
given  by  a convolution  integral.  Thus 

00 

sn(z)  - /Jj|ii(z,)P(8  - z’)dz'  (22) 

where  Sn(z)  is  the  n^  doppler  moment  of  the  scattered  pulse  shape  function 
versus  apparent  wake  station  z and  P(z)  is  the  effective  transmitted  radar  pulse 
shape.  For  backscattering  ?(z)  is  equal  to  the  relative  transmitted  pulse  power 
at  a time  (2z/c)  cos  a during  the  pulse,  normalized  such  that  the  peak  value  is 
unity.  In  this  model,  the  principal  lobe  of  a cosine  squared  function  is  used 
to  represent  E(z).  Its  width  is  defined  in  terms  of  the  parameter  range  resolu- 
tion, which  is  the  distance  between  half-power  pointB  when  a ■ o. 


III.  DETAILED  INFORMATION  FOR  THE  USER 


A.  General 

This  computer  program  was  developed  and  checked  out  using  the  Honeywell 
6060  at  the  General  Electric  Company,  Space  Systems  Division  Computer  Center, 
Valley  Forge,  Fa.  A complete  deck  of  punched  cards  (in  26-punch  and  fully 
Interpreted  In  column  alignment)  has  been  supplied  to  the  contract  monitor 
(RADC/ETEF).  A corresponding  listing  is  included  herewith  in  Appendix  A.  This 
listing  contains  the  same  statements  as  the  punched  cards.  The  deck  and  the 
listing  are  both  in  the  form  developed  for  the  Honeywell  6060,  however  no 
control  cards  are  included.  This  subsection  gives  certain  general  information 
about  the  program,  including  how  to  revise  its  central  memory  requirements  to 
do  the  problem  at  hand  economically  and  how  to  make  it  run  on  the  CDC  computers. 

The  language  is  FORTRAN-Y  and  complex,  single  precision  arithmetic  is  used. 
The  only  kind  of  type  statement  used  is  COMPLEX;  all  integer  variables  are  typed 
by  the  initial  letter,  from  1 through  N,  of  their  names.  No  mixed  mode  arith- 
metic expressions  involving  complex  variables  are  used,  despite  the  fact  that 
the  Honeywell  6060  FORTRAN-Y  compiler  permits  them,  in  order  to  avoid  any  ambi- 
guity for  the  user  and  to  avoid  difficulties  with  other  compilers.  All  state- 
ment numbers  increase  as  one  reads  down  in  each  routine,  and  none  is  super- 
fluous. Format  statements  are  all  numbered  1000  or  greater.  Output  formats' 
contain  at  least  one  space  or  a new  page  at  the  beginning  of  each  line,  to  leave 
room  for  carriage  control,  even  though  the  Honeywell  6060  does  not  need  it. 

With  certain  exceptions  noted  below  the  program  calls  only  the  built-in  and 
library  functions  iisted  in  Table  1. 

All  common  blocks  are  segregated  according  to  variable  type  and  labeled 
with  a name  which  indicates  the  typo  as  well  as  the  routines  containing  the 
common.  The  first  letter  of  the  label  indicates  the  variable  type  (C  for  com- 
plex, 1 for  integer,  and  R for  real)  with  one  exception,  ABCKLW,  which  would 
be  too  long  if  the  R were  prefixed.  The  other  letters  in  the  label  denote 
the  routines  in  which  the  common  appears.  The  code  letter  of  each  routine  ap- 
pears in  column  8 of  its  initial  comment  card.  The  13  subroutines  are  denoted 
a through,  m,  and  the  main  program  is  denoted  w for  "wake".  Thuo,  for  example, 
C0MM0N/RBCW/  consists  of  real  variables  and  appears  in  subroutines  PROFS  and 
WFP  as  well  as  in  MAIN. 

Most  of  the  dimensions  for  array  variables,  and  hence  to  a certain  extent 
the  central  memory  requirements,  are  dictated  by  the  choices  of  certain  more  or 
less  arbitrary  limits  on  the  problem  description.  The  coding  is  designed  to 
make  it  very  easy  to  change  these  limits.  The  limiting  factors  are  the  number 
of  layers,  the  number  of  angular  modes,  the  number  of  axial  stations,  and  the 
detail  with  which  radial  profiles  can  be  input.  The  radial  distribution  of  wake 
properties  is  currently  represented  by  a series  of  50  concentric  cylindrical 
layers.  The  variable  LAYERS  which  appears  in  MAIN,  FIT,  PROFS,  RINT,  WFP,  BCO, 
ST0RF,  PHINT,  FFSUM,  and  AVGS  has  this  value.  The  dimension  of  any  subscript  in. 
the  entire  program  which  now  has  the  value  50  also  is  determined  by  this  quan- 
tity. The  data  statements  for  the  arrays  ROR  and  YIP  must  of  course  conform  to 
this  number;  but  eee  the  detailed  description  below  (subsections  IIB2a)  und 
IIB2c)).  Note  that  YIP  is  dimensioned  one  less  than  this  number. 


Table  1.  Mathematical  Functions  Called  By  The  Program 


NUMBER 

TYPE 

OF 

OF 

NAME 

DEFINITION 

ARGUMENTS 

ARGUMENT 

FUNCTION 

ABS(A) 

AIMAG(C) 

ALOG(A) 

ALOGIO (A) 

CABS (C) 
CMPLX(A,B) 

COS  (A) 
CSQRT(C) 

EXP  (A) 

FLOAT  (I) 

IABS(I) 

IFIX(A) 

MAXO (l,J, . . . ) 

MOD(l,J) 

REAL(C) 

SIN (A) 

SQRT(A) 


Imaginary  part  of  C 
Natural  log  of  A 
Common  log  of  A 
|C| 

A+iB 

cos (A),  for  A In  radians 

rc~ 

0A 

Convert  Integer  to  real 

III 

Truncate  real  to 
Integer 

Maximum  value  of  all 
arguments 

I minus  the  product  of 
J with  the  truncated 
value  of  (I/J) 

Real  part  of  C 

sin (A),  for  A in  radians 

nr 


real 

complex 

real 

real 

complex 

real 

real 

complex 

real 

Integer 


real 

real 

real 

real 

real 

complex 

real 

complex 

real 

real 


Integer  Integer 

real  integer 

Integer  integer 

Integer  Integer 

complex  real 


Tha  m fields  in  the  naan  wiki  era  currently  represented  by  up  to  64 
an|ular  nodes.  Thus  tha  varlabla  MODES  haa  this  value,  which  auat  ba  an  inte- 
gral power  of  two,  in  subroutines  WFP  and  BANK.  All  variables  which  era 
how  dimensioned  64  corraapond  to  thin  requirement . , Tha  varlabla  CN(65)  in 
aub routine  HANK  auat  ba  dlnanaionad  at  laaat  one  greater  than  K0DR8.  Tha  dlnen- 
aion  2S6  for  tha  variables  0,  H and  WW  in  aubroutlnaa  FFSUM  and  AVGS  must  La 
four  time*  the  value  of  MODES.  In  AVGS  tha  vanieblala  BOS,  U and  V are  dimen- 
sioned one  greater  than  twice  MODES. 

Currently  one  can  input  up  to  20  axial  a tat ion a,  which  relates  to  all 
quantities  now  diaanaionad  20.  Note  alao  that  tha  data  statement  for  MUTZ  in 
MAIM  auat  conform  to  tha  naxiaua  number  of  axial  stations , which  is  denoted 
ISTAT  in  subroutine  21MT.  Also/  in  ZINT,  tha  variable  P nuat  ba  dlnanaionad 
at  laaat  19  or  ona  lass  than  ISTAT , whichever  is  greater.  The  quantities 
dlnanaionad  240  relate, not  only  to  tha  naxiaua  nuaber  of  axial  stations  but  also 
to  tha  fact  that,  currently,  a given  radial  profile  function  can  be  input  In 
terms  of  up  to  eleven  numbers.  Thus  the  dimensions  of  all  quantities  now  di- 
mensioned 240  must  be  the  naxiaua  nuaber  of  axial  stations  times  one  more  than 
tha  nuaber  of  quantities  Which  may  be  used  to  specify  an  input  radial  profile. 

It  it  10  daisirhd  to  input  the  radial,  profiles  in  terms  of  more  than  eleven 
numbers,  one  must  alio  change  the  dimension  of,  and  tha  data  statement  for,  the 
variable  0 In  subroutine  FIT  . (see  subsection  1112*))*. 

This  program  should  require  only  alight  modification  to  make  it  run  on 
machines  other  than  tha  Honeywell  6060.  A preliminary  version  haa  bean  convert- 
ed and  run  successfully  on  tha  CDC6600  at  APOL,  Hanscom  ATE,  MA.  Tha  forerunner 
of  this  program  waa  run  extensively  on  tha  CDC7600  at  tha  Army  BMDATC-ARC  in 
Huntsville,  AL.  The  CDG  (FORTRAN  extended)  compiler  is  a bit  more  particular 
than  the  Honeywell  (FORTRAN-Y) , in  that  certain  changes  had  to  ba  made  to  avoid 
warnings  or  fatal  error*  where  none  nod  previously  occurred.  All  auch  changes 
era  compatible  with  FORTRAN- Y and  they  have  already  bean  incorporated  into  the 
coda  as  a general  policy. 

Certain  statement » (in  MAIN,  SCO,  STORF,  end  FFSUM)  must  be  omitted  to  make! 
tha  program  run  on  tha  CDC  machines.  All  auch  statements  era  identified  in 
columns  73  to  80  by  a comment,  which  denotes  tha  reason  for  their  praaanca  in 
the  Honeywell  version.  Thie  computer  usee  a word  of  only  36  bits  limited  to 
magnitudes  between  approximately  iO*1*.  Thus  exponent  underflow  occurs  occa- 
sionally, although  it  has  s negligible  effeot  on  tha  accuracy  of  the  output. 

The  library  subprogram  FX0PT(67 ,1,1,0)  suppresses  tht  printing  of  this  srror 
massage  and  proceeds  ss  usual  with  execution.  A second  peculiarity  In  the 
Honeywell  version  is  the  use  of  scratch  random  disc  for  storage  of  certain 
quantities.  Among  tha  reasons  for  this  are  tha  facts  that  the  number  of  loca- 
tions required  depends  on  tha  input  and  that  tha  Honeywell  execution  cost  is 
proportional  to  both  oantval  memory  and  tins,  while  tha  latter  is  hardly  affect- 
ed by  disc  access,  Tha  library  subroutina  RANSIZ(lfn,n)  defines  tha  local  tile 
lfn  at  a random  disc  having  n words  par  record.  Tha  read  and  write  functions 
are  RlAD(lfn'h)liat  and  WRIT! (lfn' n)list,  where  n is  the  location  of  tha  record. 


* A caveat  is  in  order  regarding  any  changes  to  FIT  which  involve  either  tha 
"nunbar  of  inputabla  profile  points  or  the  number  LAYERS.  If  any  of  tha 
possible  input  points  is  exactly  the  same  a*  ona  of  the  value*  of  tha  inde- 
pendent variable  for  either  mash  point  or  layer  profiles,  special  programming 
techniques  nuat  ba  used  to  avoid  an  indefinite  condition. 


The  third  and  final  peculiarity  is  the  use  of  overlays  to  minimize  core  storage 
requirements*.  The  library  routine  LLlNK('name')  loads  all  routines  denoted  by 
certain  control  cards  as  belonging  to  "name"  into  an  area  where  previous  calls 
may  have  loaded  other  routines.  (The  program  hap  been  run  with  LINKA  consist- 
ing of  subroutines  WFP,  HANK,  BIN,  PROP,  BCO,  and  STORF  and  LINKB  containing 
PHINT,  FFSUM,  and  AVGS.)  In  addition  there  is  one  labeled  common  block  which 
is  common  only  to  the  two  overlaid  areas,  so  it  is  included  in  the  main  program 
to  avoid  its  loss  during  overlaying.  This  card  may  be  removed  from  MAIN  if 
overlays  are  not  used.  In  general  only  those  cards  having  an  identifier  in 
columns  73  to  80  need  to  be  removed  for  the  GDC  computers.  No  other  cards  have 
anything  in  these  columns. 

Certain  cards  must  be  added  to  make  the  program  run  on  the  CDC  machines. 

The  Honeywell  computer  does  nop  take  a PROGRAM  card.  Also,  a common  block  must 
be  defined  to  replace  the  variable  stored  on  disc.  The  following  is  a suggested 
program: 

PROGRAM  WAKE ( INPUT , OUTPUT , TAPE5- INPUT , TAPE6-0UTPUT) 

COMMON / RGH JN / E E (x ) 

CALL  MAIN 

STOP 

END 

The  quantity  x must  be  specified  according  to  the  inputs  for  the  case  being  run, 
so  this  approach  permits  one  to  control  the  amount  of  central  memory  required 
without  recompiling  the  entire  program.  The  value  of  x must  be  at  least 
12*N*(LAYERS+1) , where  N is  defined  by  the  fourth  comment  card  of  MAIN.  The 
present  main  program  should  be  headed  by  the  card  SUBROUTINE  MAIN;  and  the  card 
COMMON / RGHJN/ EE ( 1 ) should  be  placed  in  subroutines  BCO,  STORF,  and  FFSUM. 

Also,  the  appropriate  do-loop,  transferring  the  array  E into  the  array  EE  (or 
vice  versa)  replaces  the  statements  writing  (or  reading)  E on  disc.  For  example, 

L - 12*NFILE  - 11 

DO  300  1-1,12 

EE (L)  - E(I) 

300  L - L + 1 

could  be  used  in  both  BCO  and  STORF  in  place  of  the  card  identified  by  "DISC", 
in  columns  73-76  although  the  index  I need  run  only  to  8 in  BCO.  If  the  middle 
statement  is  reversed,  the  same  loop  could  be  used  in  FFSUM.  It  is  estimated 
that  approximately  100000  octal  central  memory  locations  will  be  needed  for  this 
program  in  its  present  form  on  the  CDC  machines,  not  counting  the  Bpace  needed 
for  the  EE  array.  The  latter  could  be  anywhere  from  3454  to  114400  octal  as  a 
function  of  the  inputs,  provided  that  50  layers  and  up  to  64  inodes  are  used. 


* The  program  requires  30K  memory  in  the  Honeywell  6060  when  disc  and  overlays 
are  used,  where  one  K is  defined  as  1024  decimal  locations. 
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Although  It  is  probably  not  necessary  to  do  so,  one  may  want  to  change  the 
values  of  certain  threshold  quantities  which  are  uced  to  guard  against  errors 
like  exponent  overflow  or  exponentials  with  outsit  argument,  in  view  of  the 
fact  that  10 >#  is  roughly  the  limit  for  the  Honeywell  6060.  Such  quantities  are 
THRESH  in  MAIN  and  subroutir-os  HANK,  AVGS,  FFSUM  and  ZINT;  EPSO,  STP  and  TEST  in 
subroutine  WFP}  DBT  in  subroutines  FFSUM  and  ZINT;  and  EXA  in  subroutine  AVGS. 

B.  Description  of  the  Program 

In  this  subsection  we  give  detailed  accounts  of  the  mathematical  and  logi- 
cal techniques  used.  The  listings  in  Appendix  A may  be  followed  along  with 
these  accounts . 


1.  Main  Program 

The  principal  function  of  this  routine  1b  to  direct  the  calculations  by 
calling  the  major  subroutines,  WFP,  which  calculates  the  ero  fields  in  the  mean 
waka,  PHINT  and  RINT,  which  perform  the  integrations  over  p and  <t>  in  (21),  and 
ZINT,  which  does  the  convolution  (22).  Another  function,  which  occupies  a large 
portion  of  the  executable  statements,  is  to  generate  the  radial  profiles  of 
collision  frequency,  electron  density,  and  velocity  on  the  standard  grid  of 
layers  and  boundaries  between  layers.  The  logic  is  very  simple.  There  are  two 
principal  nested  loops.  The  outer  one  is  an  implicit  loop  over  input  cases, 
between  statement  number  1 and  the  statement  immediately  before  the  STOP  state- 
ment. The  inner  loop  iterates  over  all  axial  stations.  Let  us  now  go  through 
the  listing  explaining  each  part  in  order. 


After  reading  end  writing  the  title  and  input  data,  certain  constant  quan- 
tities are  calculated.  CAY  represents  k,  the  free-space  wavenumber  (per  meter). 
WO  and  VO  represent  the  absolute  values  of  the  sine  and  cosine,  respectively,  of 
the  aspect  angle  a.  CS,  SS,  SC,  and  SMC  represent  certain  combinations  of  the 
sine  and  cosine  of  the  polarization  angle,  which  are  needed  in  subroutine  AVGS. 
EMC  is  the  critical  electron  density  Nc  (par  cubic  centimeter). 


In  the  loop  over  the  axial  stations,  EMNCO  is  the  axial  mean  electron  den- 
sity divided  by  critical,  N*/Nc,  RW  is  the  radluc  of  the  wake  r„,  CAYR  is  krw, 
and  both  VCM  and  VCOM  are  the  axial  collision  frequency  ratio  v/u. 


Shkarofsky1*  has  defined  a generalized  spectrum  function  for  isotropic  tur- 
bulence which,  when  defined  by  (20),  has  the  following  equations: 


S(k) 

U - 


(2u)*^ri<ri/r0)%(1J> 
u\(ri/*o> 
(ri/r0)/l  + (kr0)z ' 


(23) 

(24) 


In  these  equations  r-j.  and  r are  the  inner  and  outer  scales  of  turbulence,  re- 
spectively, Ky(x)  is  the  modified  Bessel  function  of  the  second  kind  of  order  y, 
and  the  parameters  y and  ip  are  interrelated  with  the  parameter  m,  which  defines 
the  shape  of  the  spectrum  function. 

if!  IP.  Shkarofsky,  "Generalized  Turbulence  Space  - Correlation  and  Wave- 
Number  Spectrum  - Function  Pairs",  Can.  J.  Phys.,  Vol.  46,  No.  19,  pp. 
2133-2153,  1 October  1968. 


.11'  - -I-int.'*..  ’-/-ill. 
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Y - (W  “ D/2  (25) 

* - (y  + 2)/2  (26) 

Now  in  order  to  obviate  both  the  need  for  calculating  modified  Beesel  functions 
and  the  evaluation  of  an  Indeterminate  form  when  the  inner  scale  vanishes,  we 
have  used  the  following  equation  in  place  of  (23): 

8ff/nr»r0|0/r(v) 


[1  + <kr0)2]'W(kr±/TT)*] 

A comparison  between  (23)  and  (27)  is  shown  in  Figure  2. 

The  program  has  the  spectrum  function  (27)  built  in  for  six  values  of  y. 

The  input  quantity  MUTZ  and  its  representative  MUT  in  the  loop  over  axial 
stations  are  defined  as  3(y  - 1).  Thus  one  can  specify  y from  4/3  to  3 in 
steps  of  1/3.  The  Kolmogorov  ppectrum  function,  which  is  obtained  by  default, 
is  represented  by  y - 5/3.  The  variable  GAMR  and  its  representative  GAR  are  the 
quantity  v9r(‘J>)/r(Y) . The  meanings  of  the  other  FORTRAN  names  for  the  spectrum 
function  parameters  are  obvious. 

One  can  sea  from  the  listing  that  the  five  sets  of  profile  functions  are 
each  generated  in  the  same  manner.  An  integer,  for  example  LC  for  collision 
frequency,  keeps  count  of  the  location  in  the  input  array  where  the  data  for  the 
profile  begin.  If  the  initial  data  point  is  zero  the  process  is  skipped,  leav- 
ing the  profiles  as  specified  at  the  previous  axial  station.  Most  of  the  pro- 
files are  generated  both  for  mesh  points  (IND>0)  and  for  layers  (IND<0),  the 
need  for  which  will  be  apparent  below. 

At  statement  number  6 the  variables  WS  and  EMNC  represent  the  values  of 
sin2ot  and  Ne/Nc  to  be  used  by  subroutine  WFP.  The  reason  for  these  new  variables 
is  that  subroutine  PROFS  changes  their  values  and  also  the  value  of  VS,  which 
is  cos2a.  TZCS  is  used  by  subroutine  AVGS. 

Statement  numbers  8 through  70  represent  output  of  the  parameters  for  the 
current  axial  station  including  all  profile  functions.  After  the  wake  field 
penetration  calculation,  MINI  (and  also  L)  Is  the  number  of  the  innermost  layer 
reached,  all  em  fields  in  layers  within  It  being  set  equal  to  zero,  AN  is 
[Na/(Ncsin2a)]2,  and  BN  is  the  denominator  of  (27)  for  S (2k) , where  k is  the  free- 
space  wavenumber.  Thus  BN  represents  the  reciprocal  of  the  non-constant  parts 
of  the  backacattering  spectrum  function  in  the  Born  approximation. 

After  integration  over  <fi  and  p (the  comment,  "integrate  over  2RDR",  refers 
to  the  latter  integral),  the  quantities,  SO,  SI,  S12,  and  S2,  represent  the 
normalised  moments  of  the  scattering  per  unit  length  for  each  of  five  polarisa- 
tion combinations.  There  are  two  kinds  of  normalisation,  which  are  removed  by 
the  "do  80".  A represents  the  Born  approximation  for  a collisionless  uniform 
wake  of  radius  ry  in  which  the  electron  density  fluctuation  equals  the  mean  elec- 
tron density.  Also,  each  normalised  spectrum  moment,  when  calculated  by  (21), 
uses  the  appropriate  electron  velocity  profile,  all  of  which  are  normalized  by 
the  axial  wake  velocity/body  velocity,  VVBZ.  The  absolute  first  moment, as  seen 
by  the  radar,  requires  multiplication  by  the  cosine  of  the  aspect  angle  ae  well 
as  by  this  normalizing  velocity.  The  second  moment  comes  from  cos2a  times  the 
"first  squared"  moment  pluB  the  true  second  moment  of  the  fluctuations . The 


quantity  S12,  representing  the  "fir At  squared"  moment,  is  the  result  of  inte- 
grating over  the  square  of  the  mean  velocity  profile.  In  the  "do  90"  loop,  the 
moments  are  converted  into  velocity  mean  and  spread,  where  the  square  of  the 
latter  is  defined  as  the  difference  between  the  mean  square  velocity  (second 
moment  divided  oy  the  zeroth  moment)  and  the  square  of  the  mean. 

2.  Subroutines 


In  this  subsection  we  follow  an  outline  in  which  the  nomenclature  is  exact- 
ly the  same  as  used  in  the  program  listing;  l.e.,  the  code  letter  for  each  sub- 
routine denotes  the  paragraph  heading  for  its  description  here. 

a)  Subroutine  FIT(N,M,IND,F,Y)  has  the  purpose  of  fitting  the  input 
data  to  a profile  on  a grid  of  either  layers  or  layer  boundaries. 

N is  the  number  of  points  input,  M is  the  location  in  the  input 
array  where  the  last  of  these  N points  is  located  (M  - N + 1 is 
the  location  of  the  firBt  data  point),  IND  is  a control  variable, 

F is  the  array  of  input  data,  and  Y is  the  array  of  output  data. 

The  array  F must  ba  dimensioned  in  the  calling  program  at  least 
sb  large  as  the  value  of  M.  IND  controls  two  aspects  of  the 
routine.  If  IND>0  the  profile  is  calculated  for  points  located 
in  the  layers;  otherwise  the  points  are  at  meBh  points  (boundar- 
ies between  layers).  If  j IND | < 2 it  is  assumed  that  the  input 
data  are  the  values  of  a smooth  function  at  N equally  spaced 
points  starting  at  the  axis  and  ending  at  the  outer  edge  of  the 
wake.  Otherwise  the  data  represent  a step  function  defined  by 
the  height  of  the  first  step  (at  the  axis),  the  relative  radius 
of  tho  first  step,  the  height  of  the  next  step,  the  relative 
radius  of  the  next  step,  etc.  Except,  of  course,  for  Y,  the 
values  of  all  variables  in  the  calling  sequence  are  returned 
unaltered . 


A second  function  of  this  routine  is  to  define  the  array  ROR, 
which  represents  the  relative  radius,  p/rw,  of  each  layer  bound- 
ary. Note  that  the  layers  are  not  uniform  in  thickness.  In 
fact  they  are  defined  so  as  to  place  a definite  ceiling  qn  the 
relative  thickness  of  Any  given  layer.  Thus,  except  of  course 
for  the  axial  layer,  the  relative  thickness,  (Pi+i~Pi)/P|  is 
not  greater  than  0.28519903.  If  the  value  of  LAYERS  were 
changed,  the  data  for  ROR  would  have  to  be  changed.  Let  c 
represent  the  limit  on  the  relative  thickness.  (See  the  dis- 
cussion of  subroutine  PROF,  where  it  is  shown  that  c should  not 
exceed  0.3.)  The  value  of  c is  minimized  if  there  is  an  in- 
teger j for  which 

(c  + 1)^  ■ (c  + l)/c 


(c  + l)1"^,  0<i<j  + 2 


Pi-1  + b 


. j<i<n 

, i - u 


isai'iiaf  *u\i, 


where 

b - rw/(n  - J + 1/c)  (30) 

and  where  n is  the  total  number  of  layers.  Of  course  j Is  the 
number  of  layers  having  uniform  relative  thickness  and  the  (J+l)th  is 
the  first  layer  to  have  the  uniform  absolute  thickness,  b,  possessed 
by  the  remaining  layers.  In  the  present  program,  j ■ 6 and  b - 
0.021049829rw.  Thus  relatively  few  layers  have  nonuniform  thickness 
and  of  the  fifty  layers  none  haB  a thickness  greater  than  about  2. IX 
of  the  wake  radius. 


Figure  3 gives  a block  diagram  of  this  routine.  As  indicated  by  the 
comment  card,  the  location  of  the  independent  variable  for  profiles 
belonging  to  points  in  the  layers  is  defined  to  be  at  the  midpoint  of 
p2  for  each  of  the  layer  boundaries  (rms  value  of  the  radius  of  the 
layer) . Thie  convention  more  or  less  guarantees  conservation  of  line 
density  in  the  layered  cylinder  approximation  of  a smooth  profile. 
Note  also  that,  for  profiles  on  mesh  points,  the  last  data  point  is 
automatically  tha  value  of  tho  last  profile  point. 
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As  tha  comment  card  indicates , N-point  Lagrange  Interpolation  is 
used  to  give  the  output  values  of  Y when  the  input  represents  a 
smooth  function  of  the  independent  variable  X.  This  approach  was 
adopted  only  after  abandoning  one  which  consisted  of  a least  squares 
fit  to  a polynomial  of  degree  (N  - 1).  The  least  squares  approach 
encountered  numerical  di . lculties,  resulting  in  rapid  divergence  of 
the  result,  when  more  than  seven  points  were  used  to  define  the  in- 
put function.  This  difficulty  occurred  even  when  the  input  function 
was  actually  an  Nth  degree  polynomial.  The  Lagrange  interpolation 
technique,  on  the  other  hand,  is  numerically  quite  stable.  For 
example  when  the  input  represents  a typical  Gaussian  at  eleven 
points,  the  output  is  no  more  than  one  unit  off  in  the  third  decimal 
place  at  any  of  the  values  of  the  array  ROR. 

The  theory  of  the  Lagrange  interpolation  for  uniformly  spaced  data 
is  right  out  of  Abramowitz  and  Stegun18.  A bit  of  manipulation  re- 
duces those  formulas  to  the  following: 

N 

Y-E  FjP(X)/(GjAj(X)]  (31) 

J - 1 


I 


'i 


ENTER 


RETURN 


FIGURE  3.  SIMPLIFIED  FLOW  DIAGRAM  FOR  SUBROUTINE  FIT. 
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The  quantity  G,  is  a symmetric  function  of  J about  the  midpoint  of 
1,. . ,,N,  where 

Gj  - (J  - 1)! (N  - J>1  , J < 1 + M/2.  (35) 

Excapt  for  P(X)  the  FORTRAN  names  follow  the  above  notation.  The 
array  G takas  advantage  of  the  symmetry  of  Gj  and  stores  the  data 
from  (35)  in  triangular  form.  The  numbers  l(2,2»3t3»A»4,5>5,6  count 
off  the  locations  of  Gj  for  N * 2,... ,11,  respectively.  The  vari- 
able L is  the  location  immediately  prior  to  the  first  one  needed 
in  a given  calculation.  The  "do  10"  loop  coverB  all  values  of  the 
independent  variable  X defined  above.  The  symmetry  of  G is  utilized 
in  breaking  up  the  sum  into  the  sequential  "do  7"  and  "do  9"  loops. 
The  variable  B is  used  to  accumulate  the  dependent  variable  Y, 
which  is  not  permitted  to  be  negative  since  it  represents  a profile 
function  for  a non-negative  physical  quantity.  The  routine  could 
in  principle  give  division  by  zero  at  the  two  locations  where  D is 
set  equal  to  F divided  by  the  quantity  G times  A.  This  will  happen 
if  any  of  the  values  of  X is  precisely  equal  to  the  value  corre- 
sponding to  one  of  the  uniformly  spaced  input  data.  This  cannot 
happen  in  the  present  program  because  the  data  for  ROR  and  the  rms 
values  of  adjacent  pairs  are  not  rational  fractions  for  numbers  up 
to  ten.  However,  it  la  not  likely  to  happen  even  if  LAYERS  were 
changed  as  long  as  equations  (28)  to  (30)  are  used  to  define  the 
values  for  ROR  and  c is  not  allowed  to  exceed  0.3.  Obviously,  of 
course,  even  if  it  happens,  one  would  simply  code  the  routine  to 
take  the  input  data  value  (Y  - F)  at  such  singular  points. 

The  step  function  routine  simply  searches  for  X values  which  fit 
on  a given  step  (the  if  statement  in  the  do-loop),  setting  Y 
equal  to  the  height  of  the  step  for  all  those  that  do,  The  steps 
can  go  up  or  down  or  both. 

b)  Subroutine  PROFS (IRE, KR1T,MHT,PSI)  has  the  primary  purpose  of 
computing  the  effects  of  scattering  by  turbulence  on  the  propa- 
gation of  the  coherent  wave  in  the  mean  wake.  Other  funotions  of 
this  routine,  which  are  included  because  of  their  logical  rela- 
tion to  either  this  function  or  the  point  in  MAIN  where  PROFS  is 
called,  are  the  coefficient  of  depolarization  by  second  scattering, 
the  adjustment  of  the  aspect  angle  and  electron  density  to  make 
the  critical  points  lie  at  layer  boundaries,  and  the  profile  of 
the  reciprocal  of  the  dielectric  constant  at  layer  boundaries • 

The  variables  IRE,  MUT  and  PSI,  in  the  calling  sequence,  are 
Input*,  which  are  unaltered  by  the  routine.  IRE  ■ 1 implies  that 
the  profile  of  mean  electron  density  la  a smooth  function  and 
that  readjuatment  of  the  axial  electron  density  and  the  sepsct 
angle  is  necassary  for  reasons  given  below.  MUT  and  PSI  pertain 
to  the  turbulence  spectrum  function  in  precisely  the  same  way 
as  in  MAIM.  KRIT  is  returned  as  zero  if  IRE  f 1 or  if  there  Is 
no  point  at  which  the  real  part  of  the  adjusted  complex  dielectric 
constant  vanishes.  Otherwise  it  is  the  number  of  the  layer 
boundary  at  which  this  critical  point  occurs.  The  logic  of  this 
routine  is  very  simple.  It  starts  with  a do  loop  over  til  layers, 
in  which  the  scattering  attenuation  and  depolarization  are 


calculated  and  the  effect  of  scattering  attenuation  on  the  mean 
plasma  Is  accounted  for.  Then  the  adjustment  of  aspect  angle 
and  electron  density  are  made,  if  required,  in  two  do  loops  over 
the  layers,  Finally  there  is  a do  loop  generating  the  profile 
of  the  reciprocal  of  the  dielectric  constant  at  layer  boundaries 
and  placing  the  adjusted  values  of  plasma  properties  in  the 
appropriate  arrays  for  use  by  subroutine  WFP.  The  five  tasks 
covered  by  the  last  three  sentences  are  each  denoted  by  comment 
cards. 

The  coefficients  of  attenuation  by  scattering  and  of  depolari- 
sation by  second  Born  are  calculated  from,  for  example,  the 
theory  given  by  Shkarofsky14.  Thus 

rft2Va  cos&n  , 

a « / (1  + x2)S(ke/2(l-x)')dx  (36) 

4[1  + (v/ui)2]  -1 

25trre2Ne' 2 1 

0 ■ /x® (1-x2 ) 2S (2kex)dx  (37) 

1 + (v/w)2  0 

where  a8  is  the  scattering  attenuation  and  0 is  the  depolarisation 
rate. *>  In  the  argument  of  the  spectrum  function  S,  the  symbol  ke 
represents  the  effective  wavenumber  in  the  mean  medium,  which  could  be 
evaluated  as  half  the  gradient  of  the  argument  of  Eo'Go.  But  this 
would  lead  to  an  iterated  problem,  so  we  let  ke  equal  the  free-space 
wavenumber  k times  the  real  part  of  the  square  root  of  the  complex 
dielectric  constant.  In  (36)  the  angle  0m  represents  a cutoff  angle 
for  omitting  the  effect  of  forward  scattering  from  the  extinction 
coefficient  as.  The  idea  for,  and  the  form  of,  0m  are  taken  from 
duthart  and  Graf15.  Thus  we  use  the  following  equation! 

COS0JJ  - 1 - [1  + (k  + ke)(rw  - p)/(2TTsina)]"1  (38) 

where  ke  is  k times  the  real  part  of  the  square  root  of  the  complex 
dielectric  constant  at  the  local  value  of  the  radius  p in  the  mean  wake. 

Of  course  the  spectrum  function  (27)  is  used  in  the  above  integrals, 
except  that  the  exponential  cutoff  involving  the  inner  scale  r^  is 
replaced  by  a sharp  cutoff.  The  integrals  reduce  to  the  following: 

X , -ill 

aB/k  - D/(16tt)/(1  - 2x  + 2x2)(l  + ax)  dx  (39) 

8 Y 

0 - kD?x2(l~x)2(l  + ax)"^dx  (40) 

0 


14 . I . P . Shkarofsky , "Modified  Born  Backscattering  from  Turbulent  Plasmas; 
Attenuation  Leading  to  Saturation  and  Cross-Polarization",  Radio  Science, 
Vol.  6,  Nos.  8,  9,  pp.  819-831,  August-Saptember  1971. 

15.  H.  Guthart  and  K.  Graf,  "Attenuation  Due  to  Scatter  in  Random  Media", 
Physics  of  Fluids,  Vol.  15,  No.  7,  pp.  1292-1299,  July  1972. 

ft  Note  that  0 is  the  two  way  depolarization;  l.e. , it  is  twice  Shkarofsky's  0. 
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»^rfr<i)>><2kr0)3Nef  2 

(41) 

r<Y)Nca[l  + (v/ui)2  ] 

(1  - cos&m)/2  , a - (2ker0)* 

(42) 

The  upper  limit  X is  equal  to  either  unity  or  the  quantity,  [it/ (2kerj) ] 2 , 
whichever  is  smaller,  the  latter  quantity  representing  the  cutoff  at 
the  inner  scale.  The  integrals  in  (39)  and  (40)  are  evaluated,  if  the 
quantity  a is  not  too  small,  by  repeated  integration  by  parts  to  the 
point  where  the  polynomial  in  x is  reduced  to  a constant  by  the  re- 
peated differentiations.  These  formulas  are  represented  in  the 
FORTRAN  statements  prior  to  statement  number  2,  and  they  differ  when 
i|i  is  integral  (MUT  - 3),  as  at  statement  number  1.  These  equations 
lead  to  indefinite  conditions  when  a is  small,  so  in  that  caiae  we 
Integrate  by  replacing  the  binomial  in  ax  by  (1  - ai|tt) , which  is  done 
between  statements  2 and  3.  In  all  statements  FJ  denotes  the  integral 
In  (39)  and  GJ  denotes  that  in  (40).  The  quantity  D at  statement 
number  3 is  D in  the  above  formulas.  AK  represents  aB/k  and  BOR, 

&/ (ksina) . 

The  attenuation  owing  to  scattering,  aB,  is  a concept  from  geometric 
optics  theories.  As  such  it  is  not  easy  to  fit  to  the  present  theory, 
which  is  a full  wave  theory.  In  optical  theories  one  integrates  as 
along  ray  paths  to  get  the  total  scattering  attenuation  to  a given 
point,  which  is  in  close  analogy  with  the  phenomenology  of  colllsional 
attenuation.  However  in  this  theory  we  get  the  local  field  at  a given 
point  by  rigorous  expansions  of  solutions  of  Maxwell's  equations.  In 
order  to  use  a representation  of  scattering  attenuation  which  is  con- 
sistent with  this  theory,  we  adopt  the  approach  of  adding  to  the  mean 
plasma  properties  in  such  a way  that  scattering  attenuation  is  auto- 
matically accounted  for.  Ideally  it  would  be  preferable  simply  to 
add  a quantity  to  the  collision  frequency  since  collisions  cause 
attenuation,  but  it  can  be  shown  that  this  is  not  possible  to  do  in 
a consistent  way.  Therefore  we  add  quantities  to  both  the  electron 
density  and  collision  frequency.  These  quantities  must  give  an 
attenuation  coefficient  which  is  Identical  with  ots  and  they  must  van- 
ish when  otg  vanishes.  There  are  many  different  sets  of  equations 
which  have  these  properties,  but  we  use  the  following: 


G - 1/(1  + ^7it) 

(43) 

D - 1 +(a8/k)2  - Ga 

(44) 

\)/w  “ 2GtxB/(kD) 

(45) 

Ne/Nc  - [1  + (v/ui)a]D 

(46) 

These  equations  are  used,  down  to  statement  number  6,  to  give  EMS  and 
CFS,  which  represent  normalised  profile  function  arrays  for  these 
added  electrons  and  collisions.  The  arrays  EM  and  CFL  are  temporarily 
used  for  the  total  profile  functions  at  mesh  points  (layer  boundaries). 

In  previous  studies  it  was  found  that  wave  propagation  in  a layered 
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plasma  can  often  be  seriously  affected  by  the  method  In  which  the 
layers  are  chosen  to  represent  a given  smooth  distribution.  This 
effect  Is  related  to  the  fact  that  the  smooth  distribution  has 
critical  density  at  only  an  isolated  point,  while  the  layered  approx- 
imation, if  selected  at  random,  might  have  critical  density  in  an 
entire  layer.  It  was  found  that  all  spurious  effects  of  layering 
could  be  avoided  if  the  two  critical  density  points  are  made  to  lie 
at  layer  boundaries.  These  two  critical  densities  are  defined  as 
that  for  which  the  real  part  of  the  dielectric  constant  K vanishes 
(intrinsically  critical  density)  and  that  for  which  the  real  part 
of  K is  equal  to  cos2a  (critical  density  for  refraction).  The  "do 
10"  loop  adjusts  EMNC  to  meet  the  first  criterion  and  the  "do  15" 
loop  adjusts  WS  (dsflned  as  slnza)  to  satisfy  the  second.  Addition- 
al conditions  satisfied  by  these  loops  are  that  the  electron  density 
be  a decreasing  function  of  radius  p at  each  critical  point  and  that 
the  outer  edge  be  underdense. 

It  will  be  seen  below  that  the  radial  component  of  the  electric  field 
is  related  to  the  tangential  magnetic  field  components  and  to  the 
reciprocal  of  the  complex  dielectric  constant.  The  arrays  ROK  and 
AOK  are  the  real  and  imaginary  parts  of  1/K  at  mesh  points . 

At  the  end  of  the  routine  the  arrays  EM  and  CFL  are  given  their 
permanent  meaning,  for  use  in  subroutine  WFP,  as  the  profile  func- 
tions for  electron  density  and  collision  frequency  in  the  layers. 

EMO  and  CFO  are  the  corresponding  quantities  without  enhancement 
owing  to  scattering  attenuation. 

c)  Subroutine  WFP(JMIN.NHP)  calculates  the  em  fields  in  the  mean  wake 
plasma*.  In  the  calling  sequence,  JMIN  is  the  number  of  the  inner- 
most layer  boundary  reached,  all  fields  within  this  point  being  assumed 
to  vanish.  NHP  is  the  total  number  of  angular  modes  UBed  to  repre- 
sent the  fields.  Both  variables  are,  of  course,  outputs  of  the 
routine. 

The  array  YIP  1b  defined  as  the  thickness  of  a given  layer  divided 
by  its  inner  radius,  starting,  of  course,  at  the  second  layer  (i.e., 
excluding  the  axial  layer) . Its  values  depend  on  the  values  of  ROR 
in  subroutine  FIT.  The  first  executable  statements  define  the  array 
RK  as  kp  for  the  outer  boundaries  of  all  layers  and  set  up  a few 
constants  for  later  use. 

The  problem  of  em  wave  propagation  in  the  cylindrical  mean  wake  plasma 
is  aolved  by  oaparation  of  variables  and  superposition.  Thus  the  vec- 
tor electric  field  E at  an  arbitrary  point  (p,<t>,a)  is  written  as  follows: 

E « axp (ikacosa)  ? e(n,p)e^n<^  (47 ) 

n«-°° 

The  quantity  e(n,p)  is  the  partial  electric  field  vector  for  the 
order  mode  and  is  a spatial  function  of  only  the  radial  coordinate. 

The  magnetic  field  H haa  the  same  representation  except  that  h replaces 


* The  name  of  this  routine  derives  from  the  name  Wake  Field  Penetration. 
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e.  In  free  space  the  functions  e and  h involve  Bessel  functions  of 
argument  kpslna  and  order  n.  The  Incident  wave  includes  only  the 
ordinary  Bessel  functions  and  the  scattered  wave  has  the  Bessel  func- 
tion of  the  second  kind  as  well.  The  asymptotic  behavior  of  these 
functions  at  large  orders  is  Buch  that  the  function  of  the  second  kind 
increases  while  the  ordinary  Bessel  function  decreases.  If  the  modal 
series  (47)  is  to  converge  at  all, its  convergence  must  depend  on  the 
ratio  of  these  two  functions.  Therefore,  in  consonance  with  the 
36-bit  computer  word  .(approximately  eight  place  accuracy)  we  truncate 
the  series  when  the  ratio  of  the  two  functions  of  argument  kr  sina 
is  greater  that  109 . In  the  section,  "determine  the  number  of  modes 
to  use",  A represents  this  argument  and  NHP  is  the  number  of  modes 
which  satisfies  this  criterion;  The  equation  for  NHP  is  a good  fit 
to  the  above  criterion  for  values  of  the  argument  between  about  0.5 
and  35,  and  It  is  conservative  at  all  other  values.  (Note  that  nega- 
tive orders  are  not  distinguished  from  positive  ones  in  this  dis- 
cussion because  of  the  symmetry  properties  which  obtain.  Thus  NHP  is 
one  greater  than  the  largest  value  of  n used  in  (47).) 

The  complex  dielectric  constant  K and  the  propagation  constant  y 


determine  the  solution  for  em  wave  propagation, 
are  used  In  the  "do  100"  loop,  are  as  follows: 

Their  forms,  which 

K - 1 - Ne/[NC(1  + iv/u))‘J 

✓-v 

00 

Y ■ k/K  - cos^ot 

(49) 

In  these  equations  Nfi  and  v are  the  electron  density  and  collision 
frequency  in  a givenlayer.  The  array  DIEL  represents  K and  the  array  \ 

GAM,  y/k.  j 

The  logic  of  this  subroutine  is  very  simple.  It  consists  of  two 
nested  loops  over  the  angular  modes  and  the  cylindrical  layers.  The 
separation  of  variables  (47)  allows  the  solutions  e and  h to  be  ob- 
tained as  functions  of  p for  each  n independently,  so  we  Iterate  over 
modes  (NP  ■ n + 1)  in  the  outer  loop  ("do  600").  In  this  loop 
EYEN=in,  SGN=(-)n,  and  ENin. 

The  solutions  for  e and  h for  a given  value  of  n come  from  substitut- 
ing both  (47)  and  its  magnetic  analog  into  Maxwell's  equations,  which 
for  exp(-iut)  time  dependence  are 

VxE  « ikH  ; V* H - - ikKE  (50) 

where  we  have  assumed  a system  of  units  in  which  the  characteristic 
impedance  of  free  space  is  unity,  for  convenience*.  These  equations 
give  six  linear  first  order  ordinary  differential  equations  in  each 
of  which  are  three  of  the  six  cylindrical  coordinates  components  of 
e and  h.  These  six  equations  reduce  to  four  by  elimination  of  the 
radial  components,  but  each  of  these  four  equations  contains  in  gene- 
ral three  of  the  four  tangential  components.  The  easiest  way  to 


* The  system  of  units  is  immaterial  since  all  results  will  eventually  be 
expressible  in  nondimensional  terms. 


25 


solve  these  ie  to  suppose  that  all  fields  are  the  sums  of  two  kinds, 

TE  and  TM,  where  T denotes  transverse  to  the  axis  and  E and  M refer 
to  electric  and  magnetic.  Thus,  for  TE,  when  ez  ■ 0 is  substituted 
one  gats,  by  elimination  of  the  <j>  - components,  a second  order  differ- 
ential equation  in  hz.  A similar  equation  for  ez  results  from 
assuming  hz  ■ 0.  These  two  equations  differ  only  in  terms  of  the 
effect  of  gradients  of  K,  so  in  our  layered  approximation  they  are 
identically  the  same  equation,  which  is  Bessel's  equation  for  order 
n and  argument  £,  where 


5SYP  (51) 

The  four  simultaneous  differential  equations  can  also  be  solved  for 
the  4>  - components  as  functions  of  the  z - components.  The  results 
are  expressible  in  terms  of  the  following  variables: 

pHik/y  (52) 

q=nkcosa/(Y£)  (53) 

Thus 


1 1 
j . 

M 


■■ 


e«j)  ■ - qez  - phz'  (54) 

h<j,  « - qhz  + Kpez'  (55) 

where  prime  denotes  differentiation  with  respect  to  £. 

Each  tangential  field  component  must  be  continuous  across  boundaries 
between  layers.  A matrix  representation  is  a convenient  way  to  keep 
satisfying  this  condition  from  layer  to  layer.  Thus  we  arbitrarily 
adopt  the  following  tangential  fields  vector  F: 


F= 


h 

h- 


(56) 


- e 


♦J 


Now,  for  the  time  being,  we  restrict  the  discussion  to  the  behavior 
of  F within  any  given  layer.  Since  the  z - components  are  each 
Independent  linear  combinations  of  two  Independent  solutions  of 
Bessel's  oquatlon,  let  F - UA,  where  A is  a column  vector  of  undeter- 
mined coefficients  and  the  square  matrix  U comes  from  the  equations 
for  each  component  of  F.  Let 
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where  u and  v are  the  twd  solutions  of  Besael's  equation.  Let  sub- 
script 0 denote  the  value  of  a quantity  at  the  inner  boundary  of  a 
given  layer,  unuubs  •ripted  variables  denoting  either  constants  or 
the  values  of  the.  \ «r;>bles  at  the  outer  edge  of  the  layer.  There 
, -art  exist  a square  matrix  P such  that  F ■ PF0,  which  implies  that 
tlA  • PUoA  since  A is  constant  throughout  the  layer.  Thus  U » PU0  or 


P * UUi;”1.  Now  if  we  require  0/  the  two  solutions  that  v0  ■-  uj  ■ 0, 


then  U0  has  only  two  nonsaro  off-diagonal  elements  and  is  easy  to 
Invert.  The  result  is 


P« 


u/ur 


V/ (Kpv0 ' ) 


[Kpu ' / u0+qq0v/ (pv0 ' ) ] v ' / v0 ' 


-q0v/(pv0') 


q0v/(Kpv0’) 

<qov,/vo'“^uo) 


-q  V/  (pv0' )] 
v/(pv  ') 


(qu/u0-q0v' /v0' ) 


qv/(Kpv0')  [qq0v/(Kpv0')+pu*/u0]  v’/v,,’ 


(58) 


Three  of  the  fourteen  nonsero  elements  of  this  matrix  are  repeated, 
so  the  arrays  P P,  PQ,  . ..,  PZ  represent  the  eleven  unique  elements  In 
their  order  of  appearance  reading  across  rowe.  The  variables  P,  Q,  R, 
and  S,  which  ara  calculated  in  subroutine  PROP,  represent  u/u0  , 


v/v( 


1 

0 * 


/ uc,  and  v'/v0't  raspactivaly. 


'0  > 


Now  the  boundary  condition  between  layers  Is  automatically  eatiafiad 
if  F la  matched.  Therefore  F at  the  outer  edge  of  a given  layer  la 
identical  with  F at  the  inner  edge  of  the  next  layer.  Thus,  since  P 
represents  the  propagation  of  F through  a layer,  cumulative  matrix 
multiplication  of  the  P matrices  for  each  layer  glvos  the  propagation 
of  F through  all  layers.  This  cumulative  product  has  ths  FORTRAN  name 
QQ(i,J)  for  the  alamant  of  the  1th  row  end  jth  column.  Initially  . 
it  la  the  identity  matrix  (the  "do  200"  loop)  end  it  is  updated  in  '.he 
"do  300"  loop.  One  will  note  that  QQ  is  replaced  by  QQ  times  P as  the 
layers  art  iterated  from  the  outside  in,  Thus,  sines  P raprastnts 
propagation  of  F from  the  inside  to  the  outside  of  the  layer,  QQ  is 
the  propagation  of  F from  the  currently  reached  inner  boundary  to  the 
outermost  boundary. 


The  overall  propagation  matrix  (1st  us  call  It  Q instead  of  QQ  In 
this  discussion)  tends  to  diverge  in  cases  in  which  there  is  strong 
cutoff  of  ths  am  fields  in  ths  csntrel  regions  of  ths  layered  cylinder. 
Owing  to  the  nature  of  the  Wronskian,  the  inverse  representation  Q“J , 
which  gives  propagation  inward  from  tha  outside  edge,  also  diverges. 
Thua,  no  matter  how  one  calculates  propagation,  there  can  be  a point 
at  which  exponential  overflow  will  occur,  especially  in  case  of  a 
limited  computer  word.  In  the  computer  modal  we  insert  a perfectly 
conducting  boundary  just  before  that  point  is  reached.  Three  criteria 
ara  used.  Tha  first  is  the  if  statement  near  the  beginning  of  the 
"do  500"  loop,  which  guards  against  tha  fact  that  the  layer  propaga- 
tion matrix  P tends  to  Increase  exponentially  with  tha  argument  £. 

Tha  second  la  the  two  if  statements  in  the  "do  400"  loop,  which  guard 
against  overflow  of  |q|2.  The  third  criterion,  which  is  exercised 
just  after  the  latter  loop,  is  based  on  the  supposition  that  t (celled 
EP80)  represents  the  smallest  field  intensity  of  any  practical 
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significance.  Then  the  greatest  value  of  Q of  any  interest  ia  such 
that  |q| 2c  has  about  the  magnitude  of  the  field  Intensity  at  the 
outer  udge,  which  is  approximately  proportional  to  the  square  of 
the  ordinary  Bossel  function  of  argument  krwsina.  The  Bessel  func- 
tion is  in  the  neighborhood  of  unity  for  small  orders  and  at  large 
orders  it  drops  off  according  to  well  known  asymptotic  forms. 

Combining  these  limits  into  a single  expression  gives  the  criterion 

|q| 2e{l  + /2fln[2n/(ekrwsinct)]n}*<l  (59) 

where  e is  the  base  of  natural  logarithms.  The  magnitude  squared  of 
Q is  estimated  from  its  average  value  for  all  sixteen  elements. 

The  section  of  this  routine  which  involves  boundary  conditions  is 
closely  related  to  subroutines  BIN  and  BCO.  The  theory  is  detailed 
under  those  routines  rather  then  here  because  this  routine  merely 
celculatee  certain  parameters  for,  or  from,  those  routines. 

d)  Subrountlne  HANK(NP,X,HNK,HNKP)  calculates  the  Hankel  functions  of 
the  first  kind  and  their  derivatives  for  use  in  the  outer  boundary 
condition.  NP  ie  the  number  of  orders  to  be  calculated  (n  - o to 
NP  - 1)  and  X is  the  argument,  a real  variable.  The  complex  arrays 
HNK  and  HNKP  contain  the  output  Hankel  functions  and  their  deriva- 
tives and  must  be  dimensioned  at  least  NP  by  the  calling  program. 

The  method  used  is  downward  recurrence  on  the  Bessel  function, 
normalization  baaed  on  the  value  unity  for  the  generating  function, 
the  Neumann  expansion  of  Y0(x)  in  terms  of  Jn(x)  for  all  even  values 
of  n,  the  use  of  the  Wronokian  to  get  Y^x),  and  finally  upward  re- 
currence on  the  function  of  the  second  kind,  Yn(x).  Refer  to 
chapter  9 of  Abramowitz  and  Stegun  (op.  cit.,  ruf.  13)  for  the 
equations,  which  are  readily  located  in  the  listing  with  the  aid  of 
the  comment  cards. 

The  form  of  the  convergence  criterion,  which  identifies  the  order  at 
which  to  begin  the  downward  recurrence,  is  based  on  the  behavior  of 
Bessel  functions,  for  which  it  can  be  shown  from  tabulated  values 
and  asymptotic  forms  that  J-(x)<10~y  if  n is  approximately  equal  to 
or  greater  than  [y  - 4.5  + /3|x| ( |x|  + 8) ] . The  value  12  fr-  y was 
derived  from  experimental  (checkout)  rune  using  the  36-bit  computer 
word. 

e)  Subroutine  BIN(X,NP,B1)  calculates  Jn'(X)/Jn(X)  for  complex  X and  for 
orders  n ■ 0 to  NP  1.  The  results  are  placed  in  the  complex  array 
Bl,  which  must  be  dimensioned  at  least  N?  in  the  calling  routine. 

The  method  used  is  downward  recurrence;  l.e., 

VOO/JnW  " n/X  - cn+l<x>/cnW  (60) 

Cn(X)  - 2(n  + l)Cn+10O/X  - Cn+2(X)  (61) 

where  the  Cn's  are  started  with  arbitrary  values  at  the  order  given 
by  the  convergence  criterion,  which  was  derived  as  in  oubroutine 
HANK.  Of  course  normalization  using  a generating  function  is  not 
required.  Also,  there  is  no  danger  of  division  by  zero  as  long  as 
the  Imaginary  part  of  X does  not  vanish  (collision  frequency  not 
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identically  taro  In  tha  axial  layer),  bacauaa  all  zeros  of  Jn(X) 
are  real. 


I 


Tha  raiaon  d'etre  of  this  routine  is  the  inner  boundary  condition. 
In  case,  In  subroutine  WFP,  the  axial  layer  is  reached  in  the 
propagation  calculation,  the  boundary  condition  at  its  surface  must 
account  for  the  presence  of  only  the  Bessel  function  of  the  first 
kind.  We  state  this  boundary  condition  in  the  following  fora  in 
terms  of  the  two  unknown  constants  Aj  and  A2 : 


Fi  - 


*1  Bn 


A*  - q At 
A* 


AaBia  + dAiBn 


(62) 


B12  * ikB^y 
Bn  - If  (B|2K) 


(63) 

(64) 


Jn’(0/J„(0 


(65) 


where  the  dielectric  constant  X and  the  quantities  y,  $ and  q,  which 
are  defined  by  (49),  (51)  and  (53),  are  for  the  axial  layer  at  its 
outer  edge.  Of  course  F^  represents  tha  tangential  fields  vector  F 
at  this  boundary.  This  normalisation  is  used,  rather  than  one  in 


which  Jn  and  Jn'  appear  separately  in  all  terms  of  (62),  because  Bi 


is  easier  to  calculate  and  because  this  fora  is  compatible  with  the 
perfectly  conducting  boundary  condition,  in  which  Bu  and  Bla  are 
leero  to  account  for  the  vanishing  of  tha  tangential  electric  fields. 


f) 


Subroutine  PROP(EN,XZR,XZI,Y,C,D,CP,DP)  calculates  the  four  basic 
quantities  which  comprise  the  nonasro  elements  of  the  propagation 
matrix  P for  a given  layer.  Tne  first  four  vsrisbles  in  the  calling 
sequence  represent  input'*  and  tha  last  four,  outputs.  EN  is  the 
floating  point  value  o c the  mods  order  n.  XZR  and  XZ1  are  the  real 
and  imaginary  parts  of  the  quantity  which  represents  yp  at  the 
inner  boundary  of  the  layer.  Y is  the  relative  layer  thickness.  For 


the  purposes  of  the  derivations,  1st  C,  D,  C , and  D'  represent 


C,  D,  CP,  and  DP,  which  are  called  P,  Q,  R,  and  S in  subroutine  WFP. 


In  the  limit  of  ssro  layer  thickness  the  propagation  matrix  must 
approach  the  identity  matrix.  In  this  limit  the  quantities  C and  D', 
which  are  on  the  diagonal,  must  go  to  unity  and  the  quantities  D and 
C'  must  vanish.  Supposing  that  this  limit  is  approached  continuously, 
wa  use  series  expansions  in  the  layer  thickness  to  calculats  those 
quantities.  The  limiting  forms  serve  to  determine  the  forms  of 
these  eeries  up  to  the  second  order.  Thus  let 
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D - $0Y[1  + Z d,] 
j-1  J 

00 

C'  - E (j  + l)Ci 
j-1  J 

00 

D'  -1+2  (j  + l)di 

j-1  J 


(68) 


(69) 


In  these  equations  cj  and  dj  are  each  equal  to  a complex  number,  which 
depends  on  j,  tlmeB  the  j ^ power  of  the  relative  layer  thickness  Y. 
The  recurrence  relation  for  these  terms  Is  derived  by  substitution  of 
either  the  C or  the  D aeries  In  Bessel's  equation » for  which  the  inde- 
pendent variable  is  £ - £0(1  + Y),  and  requiring  that  the  result  apply 
Independently  for  each  power  of  Y.  Thus 


j(j  + Daj  + (2J  - DjYaj.j  + l(j  - 1>2  - n2]Y2aj_2 
+ (aj-a  + 2Yaj_ j + Y2aj_4)  - 0 


(70) 


where  aj  represents  either  cj  or  d^ . Thla  recurrence  relation  applies 
liberally  as  long  as  the  followingJ initial  conditions  are  usad: 


■j  ■ 0 , j<  - 1 

(71) 

c0  - d«j  - 0 

(72) 

■ 1/«0Y> 

(73) 

d0  - 1 

(74) 

In  the  program  the  approximate  convergence  criterion  represents  a fit 
to  the  minimum  number  of  terms  needed  in  the  Infinite  series  to 
guarantee  convergence  in  all  of  more  than  500  runs  of  this  routine 
which  were  made  during  checkout.  These  sumB  also  have  an  Internal 
convergence  criterion  following  the  "do  8"  loop.  The  array  AJ(I,J,K) 
represents  cj  for  I ■ 1 and  dj  for  1-2.  The  real  and  imaginary 
parts  are  for  J - 1 or  2,  respectively.  The  five  locations  in  the 
third  subscript  represent  the  values  needed  for  recurrence.  The  values 
of  these  locations  are  denoted  L,  LA,  LB,  LC,  and  LD  and  are  updated 
in  each  iteration.  The  first  three  termB  of  (70)  have  real  coeffi- 
cients and  these  are  represented  temporarily  in  terms  of  the  array  P. 
The  coefficient  of  the  complex  variable  is  temporarily  placed  In 
the  array  Q.  After  updating  the  aj  the  sumB  used  in  C and  D are 
accumulated  in  U and  those  used  in  C'  and  D',  in  W,  but  not  before  P 
and  Q are  temporarily  used  for  the  old  values  of  U and  W to  facili- 
tate convergence  testing.  At  the  end  the  outputs  are  calculated  as 
in  (66)  to  (69). 


j?  The  present  form  of  this  routine  is  believed  to  be  the  optimum  for  a 

J:  multi-layered  medium.  These  four  matrix  elements  are,  of  course, 

« identical  with  the  cross  products  of  the  first  and  second  kinds  of 

!•  I Bessel  functions  for  the  complex  arguments  for  either  boundary  of  the 

? ' layer.  But  the  method  of  calculation  based  on  this  fact  would  re- 


layers.  The  recurrence  relations  for  these  cross  products  were  tried, 
but  no  way  was  found  to  eliminate  the  fatal  round  off  error  which  was 
accumulated.  A detailed  study  of  series  expansions  whose  leading 
terms  are  designed  to  fit  the  various  asymptotic  forms  of  Bessel  func- 
tions was  done  also,  however  all  of  these  failed  in  some  regime  where 
the  above  scheme  worked.  With  the  constraint  that  Y be  not  greater 
than  0.3  this  algorithm  works  in  the  Honeywell  6060  for  values  of  the 
outer-edge  argument  |$9(1  + Y)|  as  large  as  88.  This  ia  exceptionally 
good  performance  because  the  result  tends  to  be  exponential  in  this 
argument  and  exp(88)  is  at  the  limit  of  the  computer  word.  Note 
that  the  implicit  series  in  Y,  in  which  Y^  is  part  of  c*  and  dj , is 
largely  responsible  for  this  good  performance.  When  explicit  power 
series  wars  used  we  could  not  get  past  about  50,  instead  of  88,  with- 
out exponent  overflow. 

g)  Subroutine  BCO(M)  applies  the  boundary  conditions  co  get  the  local 
fields  et  all  mesh  points  from  layer  M out  to  the  outer  edge. 


The  outer  boundary  condition  is  the  sum  of  the  tangential  fields  for 
the  incident  and  scattered  waves.  Standard  Fourier-Bessel  expansions 
of  the  vector  incident  fields  given  by  (16),  (17)  and  (18)  produce 
expressions  for  the  tangential  components  of  e in  the  Incident  wave. 
The  incident  waves  are  so  defined  that  the  magnetic  field  for  trans- 
verse .polarisation  is  identical  with  the  electric  field  for  parallel 
while  the  magnetic  field  for  parallel  is  the  negative  of  the  electric 
field  for  transverse.  The  result  is 


(75) 
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in  Jn(krw  sin  o)  sin  a 


Fjj  - in  + 1 jn'  (krw  ain  a) 


(76) 

(77) 


where  F*  denotes  the  contribution  of  the  incident  wave  to  the  tangen- 
tial fields  vector,  defined  by  (56),  at  the  outer  edge  and  b is  q, 
defined  by  (53), evaluated  in  free  space  at  the  outer  edge. 


b ~ n cos  a/ (krw  *in2a) 


(78) 


The  axial  components  of  the  scattered  waves  must  be  proportional  to 
the  Henkel  function  of  the  first  kind,  so  we  let  the  scattered  tangen- 
tial fields  vector  F8  have  the  following  form: 
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Af  B2  - b A* 


Af  B2  + b A® 


Ba  - i H^(krw  sin  a)/[  H^kr*,  sin  a)  sin  a]  (80) 

where  A®  end  A®  are  unknown  scattering  coefficients. 

The  boundary  conditions  are  expressible  as  the  following  matrix  equa- 
tion: 

pi  + F°  - Q Ft  (81) 

where  one  will  recall  that  1b  the  inner  boundary  condition  given 
by  (62),  and  that  Q represencs  propagation  of  F from  the  inner  boundary 
point  to  the  outer  edge.  Because  it  involves  the  four  unknown  co- 
efficients in  F£  and  Fs,  this  equation  is  equivalent  to  another  matrix 
equation  as  follows: 


2 8 
11  12 


S21  S22 
S,l  S32 


S4!  S42 


where 


sji  ■ Bu  <Qj»  + q Qj4>  + Qj2 


Sj2  " 3 + B12  ^j4  “ 4 ^j2 


It  is  not  necessary  to  solve  this  equation  for  A®  and  A®  because  At  and 
A2  determine  the  fields  at  all  points  via  the  propagation  matrix  P for 
each  layer,  Equation  (82)  represents  the  equation  SA  » pi,  which  is 
transformed  by  the  matrix 

B2  -1  -b  0 

T - (85) 

b 0 B2  -1 

into  the  following  equation  of  rank  two: 
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The  quantity  F3  is  (B2Ft  - F2),  but  the  Wronskian  gives 


(86) 


F9  - - 2in/[7f  krw  Hn(krw  sin  a)  sin  a]  (87) 

The  two  systems  of  equations  represented  by  (86)  are  readily  solved 
for  the  two  alternative  sets  of  values  of  Aj  and  A2.  Note  that  all 
results  are  proportional  to  the  quantity  F3.  The  complex  exponential 
ein$  is  used  in  the  modal  expansion  (47)  because  it  represents  all 
field  components  simultaneously.  However  each  cylindrical  coordinate 
component  of  each  (electric  and  magnetic)  field  for  each  incident 
wave  -polarization  has  either  even  or  odd  symmetry  in  n and  Thus 
(47)  is  actually  a summation  over  positive  n involving  either  2 cos(n<}>) 
or  21  ain(n(j>)*  as  well  as  the  term  for  n - 0,  Therefore  we  take 
care  of  the  factor  2 for  n>0  by  multiplying  F3  by  two  in  subroutine 
WFF,  where  this  quantity  is  calculated.  (The  quantity  Ft,  which  is 
needed  to  give  the  coherent  scattering  coefficients,  is  likewise  doubled 
at  the  same  point.) 

Now,  to  Interpret  the  FORTRAN,  the  variable  07  represents  b and  B 
represents  q.  C5  and  C6  are  temporarily  used  for  SX1  and  S3l,  and  Cl 
and  C2  are  the  elements  in  the  firQt  column  of  TS.  Than  C5  and  C6  are 
S.2  and  S32,  and  C3  and  C4  are  the  second  column  of  TS.  The  two  sets 
of  simultaneous  equations  (86)  are  then  Bolved  and,  after  the  "ao  10" 
loop  the  arrays  FJA(K) , FJB(K),  FJC(K),  and  FJD(K)  are  the  four  ele- 
ments, in  order,  of  the  tangential  fields  vector  F.  The  first  set 
(K  ■ 1)  is  for  a parallel  polarization  incident  wave  and  the  second 
(K  ■ 2)  is  for  transverse.  After  this  the  "do  200"  loop  generates  the 
tangential  fields  at  each  layer  boundary  by  Iterated  matrix  multipli- 
cation by  the  stored  propagation  matrix  P for  each  layer. 

The  final  function  of  this  routine  is  to  generate,  output,  and  store 
(on  disc)  the  coherent  scattering  coefficients,  starting  at  statement 
number  210.  Note  from  the  definition  (56)  of  F that  according  to  (79) 
and  (75)  the  scattering  coefficients  are  fully  characterized  by  the 
z-components  of  the  fields  at  the  outer  boundary  after  Ft  has  been 
subtracted  from  ez  for  parallel  polarization  and  hz  for  transverse. 

Now  in  free  apace  the  four  z-ccmponents  each  propagate  as  Hn(kp  sin  a), 
so  their  values  at  any  point  may  be  gotten  by  multiplying  the  ratio 
of  this  quantity  to  its  value  at  the  wake  edge  by  their  values  at  the 
wake  edge.  In  the  far  field  this  ratio  has  the  asymptotic,  form 

2 exp[i(*Tr/4  -*•  kp  i»in  a)] 


/ 2nkp  sin  a in  Hn(krw  sin  a) 

n 2 

Now  see  subroutine  WFP,  where  Z is  0.5  insina  Hn(krw  sin  a).  Also 
noting  that  the  slant  range  r from  the  origin  of  the  wake  is  equal  to 
p/sin  a,  this  implies,  for  example,  that  in  the  far  scattered  field  for 
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parallel  polarization 
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j' 
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aln  a exp[i(kr  - ir/4)]  „ 

Ez  — £ C3  cob  n (J)  (88) 

/"2t r kr  n-o 

where  C3  is  the  FORTRAN  variable  which  is  output  and  stored  (on  disc) 
at  this  point  in  the  program.  Now,  seeing  that  the  scattered  wave 
propagates  in  the  direction  making  angle  ot  to  the  z-axis, 

i(kr  - tt/4)  „ 

‘‘scattered  parallel  . e _ j C3  cos  n tj)  (89) 

/ 2iTkr  n"° 

This  is  a component  of  the  coherent  far  electric  field  scattered  by  a 
wave  incident  in  the  plane  of  the  axis  (as  defined  above).  This  com- 
ponent also  lies  in  the  plane  of  the  axis,  where  the  positive  direction 
is  as  defined  for  the  incident  wave.  The  same  component  of  the  scat- 
tered magnetic  field  for  transverse  incident  wave  polarization  has 
identically  the  same  equation  with  C6  replacing  C3.  These  two  field 
components  have  even  symmetry.  The  other  two,  having  odd  symmetry, 
use  1 sin  nij>  in  place  of  cos  nij>.  The  axial  plane  magnetic  field  for 
parallel  incident  wave  polarization  has  C4  in  the  summation  and  the 
axial  plane  electric  field  for  transverse  polarization  has  C5.  The 
fact  that  half  of  these  results  are  magnetic  fields  presents  nn  prob- 
lem conceptually  since  in  the  far  field  the  electric  vector  is  equal 
to  the  cross  product  Hx$,  where  f is  a unit  vector  in  the  direction  of 
scattering.  Thus  the  transverse  (<{>  - component)  electric  field 
scattered  by  a parallel  incident  field  is  the  sum  in  C4  and  the  trans- 
verse electric  field  for  transverse  incident  field  uses  C6.  The  write 
statement,  then,  gives  the  coefficients  for  parallel  and  transverse 
in  order,  first  for  parallel  incident  field  and  then  for  transverse. 

For  reasons  of  convenience  below,  the  real  and  imaginary  parts  of  the 
field  components  are  stored  in  odd,  even  order,  where  the  odd  compon- 
ents have  the  coefficient  multiplied  by  i so  that  the  Fourier  sums  are 
over  either  sin  n<j>  or  cos  n<|>.  This  order  is  Etp,  Ett>  Ept»  Epp>  where 
the  first  subscript  is  the  incident  wave  polarization  and  the  second  is 
the  scattered.  The  storage  location  Implies  that  these  are  scattered 
fields  in  that  one  extra  "layer"  is  saved  for  this  (LAYERS  + 1). 

h)  Subroutine  STORF(NP, J,ENKR,V)  stores  in  odd,  even  order  the  local 

electric  fields  components.  All  elements  of  the  calling  sequence  are 
inputs,  which  are  unaltered  by  the  routine.  NP  is  (n  + 1),  J is  the 
number  of  the  layer  boundary,  ENKR  is  n/(kp)  for  this  boundary,  and  V 
is  cos  a.  All  components  having  odd  symmetry  are  multiplied  by  i in 
order  that  all  total  fields  be  straight  Fourier  sums.  The  order  of 
storage  is  real,  imaginary;  p,  cj> , z - component;  transverse,  parallel 
incident  wave;  where  the  first  set  of  things  varies  the  most  rapidly 
and  the  last  varies  the  least  rapidly. 

The  radial  electric  fields  are  given  by  the  following  equation,  which 
is  one  of  the  six  resulting  from  substitution  of  the  modal  expansions 
(47)  in  Maxwell's  equations: 
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e p “ { h(f,  cos  a - [n/(kp)]  hz)/K 


The  general  storage  location  denotes  both  the  layer  boundary  and  the 
mode  order, 


Subroutine  PRINT (LS, MM)  directs  the  analysis  of  the  $ - dependence  of 
all  field  quantities  for  the  mean  wake.  One  of  the  most  important 
aspects  of  these  calculations  is  the  integration  over  <}>  for  equation 
(21),  whence  the  name  PRINT.  In  the  calling  sequence.  LS  is  the  number 
of  the  innermost  layer  boundary  reached  in  the  field  penetration  cal- 
culation. MM  is  the  total  number  of  modes  used  in' the  expansions  of 
the  fields. 

The  first  part  of  this  routine  sets  up  the  grid  of  angles  to  use  for 
representing  the  variation  of  all  fields  with  The  criterion  1b 
that  the  number  of  uniformly  spaced  angleB  be  a power  of  two  which  is 
at  least  four  times  the  number  of  modeB  UBed  in  the  field  expansions. 
This  criterion  derives  from  the  fact  that  the  integral  over  all  angles 
<t>  of  any  m**1  order  Fourier  sum  is  precisely  represented  by  the  trape- 
zoidal rule  for  evenly  spaced  valueB  of  (j)  as  long  as  the  number  of 
angles  is  greater  than  or  equal  to  m.  The  product  |e  *G0|?  represents 
a Fourier  sum  of  order  4m  if  E0  and  G0  are  each  of  order  m.  In  the 
routine,  IM  is  the  number  of  angles  and  M is  the  power  of  2 to  which 
it  is  equal. 

The  routine  calls  for  the  summation  of  the  coherent  scattered  fields 
before  going  to  the  loop  over  layer  boundaries,  at  each  of  which  the 
- dependence  is  analyzed. 

Subroutine  FFSUM(LYR,MM,N,M,LTH)  UBes  the  FFT  algorithm  to  evaluate  the 
various  Fourier  sums  which  represent  the  components  of  the  electric 
field  in  the  mean  wake.  LYR  is  the  number  of  the  layer  boundary  or, 
if  greater  than  LAYERS,  an  Indication  that  the  coherent  scattered 
fieldB  are  to  be  evaluated.  MM  is  the  number  of  modes  for  which  the 
Fourier  coefficients  of  the  fieldB  have  been  calculated.  H is  the 
number  of  angles,  which  must  be  the  M*"  power  of  2,  If  LTH  - 0 the 
routine  assumes  that  sin  <J)  and  cos  $ for  all  (j)  on  the  grid  of  angles 
have  not  been  calculated  before.  These  quantities  are  placed  in  the 
array  WW  during  execution  and  LTH  • 1 is  returned.  If  LTH  ■ 1 it  is 
assumed  that  these  values  already  are  stored  in  WW.  (These  assumptions 
are  valid  as  long  as  the  calling  routine  (PHINT)  does  not  change  the 
values  of  N and  M without  resetting  LTH  * 0 and  as  long  as  the  array 
WW  is  not  changed  by  the  routine  (AVGS)  with  which  it  has  this  variable 
in  common . ) 

This  routine  is  based  in  greater  part  on  the  subroutine  FFCPLX  by 
Martin!6  In  the  present  notation,  the  arrays  G and  H represent  both 
input  and  output.  Thus  the  FFT,  between  "*bit  inversion"  and  statement 
number  40,  represents  the  following: 


16.  M.A.  Martin,  "Frequency  Analysis  of  Real  Sampled  Data  by  the  Fast  Fourier 
Transform  Technique",  General  Electric  Co.,  Technical  Information  Series, 
Report  No.  69SD267,  May  1969. 
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G(I,K)  + lH(IfK)  - Z [G(I, J)  + iH(I,J)]e1(J  " (91) 

J-l 

where  4>  ■ 2tt(K  1)/N  and  K also  runs  from  1 through  N.  The  bit 
inversion  part  of  the  routine,  in  which  the  locations  of  all  variables 
are  changed  to  the  location  whose  binary  digits  are  the  reverse,  per- 
mits the  output  to  occupy  the  same  places  as  the  input.  The  econo- 
mized Fourier  summation  is  essentially  the  Cooley  and  Tukey17  method. 

Note  that  the  imaginary  part  of  the  right  hand  side  of  (91)  is  a 
cosine  series  if  the  G's  are  all  zero  or  a sine  series  If  the  H's  are 
zero.  Thus  we  apply  the  fast  Fourier  sum  by  inputting  ell  even  func- 
tions to  the  H array  and  all  odd  ones  to  G;  the  output  is  then  in  H 
in  all  cases.  All  members  of  both  arrays  are  Initially  set  to  zero 
and  then  tho  data  for  the  modes  calculated  by  WFP  are  read  in.  The 
data  are  missing,  of  course , if  the  innermost  layer  boundary  reached 
by  WFP,  which  is  denoted  INLAY,  is  outside  of  the  point  currently 
being  considered.  The  arrays  G and  H are  doubly  subscripted  to  permit 
the  fast  Fourier  sums  to  be  evaluated  simultaneously  for  all  8 or  12 
field  components  which  were  stored  by  BCO  or  STORF.  Alternating 
pairs  of  these  components  are  odd  and  even,  ho  alternating  pairs  are 
input  to  G and  H in  the  "do  2"  loops. 

If  the  point  in  question  is  within  the  wake  boundaries,  the  routine 
proceeds  to  the  evaluation  of  the  Integral  over  <{>  Involving  IEq’GJ2 
Otherwise  it  outputs  the  coherent  scattering  versun  4),  which  goes  from 
0 to  180°,  the  other  half  of  the  space  being  symmetrical  with  this 
half.  Referring  back  to  g)  Subroutine  BCO,  the  complex  parallel 
polarization  scattering  1a  in  the  arrays  H(7,  J)  and  H(8,  J).  Like- 
wise elements  3 and  4 are' for  transverse.  Linear  orthogonal,  which 
is  X(3),  is  taken  as  the  rms  value  of  ELt  and  Etp.  Circular  polariza- 
tion is  straightforward  also  Bince  the  transmitted  field  in  composed 
of  one  linear  polarization  pluB  i multiplied  by  the  other.  We  then 
define  principal  circular  polarization  scattering  as  the  conjugate 
coubination;  l.e., 

principal  * £ l*tt  + 1 «pt  - 1 <*tp  + 1 V1  (92) 

and  orthogonal  is  therefore 

^orthogonal  “ 7 [®tt  + * %>t  + * (^tp  + i Epp)]  (93) 

These  equations  are  used  to  give  X(4)  and  X(5)  in  the  program.  Of 
some  Interest  in  calculating  far  field  diffraction  Into  the  shadow  of 
the  wake  is  the  coherent  forward  scattering.  The  complex  coefficients 
of  Epp  and  Ett  for  <j>  ■ 0 are  output,  after  accounting  for  in  (89), 

in  tne  variables  A,  B,  C,  and  D. 
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17.  J.W.  Cooley  and  J.W.  Tukey,  "An  Algorithm  for  the  Machine  Calculation 
of  Complex  Fourier  Series",  Mathematics  of  Computation,  Vol.  19,  pp  297- 
301,  1965. 
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Subroutine  AVGS(H,IK,N,LAYR)  has  the  main  purpose  of  performing  the 
Integral  in  (21)  over  all  angles  4> . The  variable  H in  the  calling 
sequence  represents  the  IK  components  of  the  electric  field  at  N 
equally  Bpaced  values  of  <j>  from  0 to  360° . IK  must  be  equal  to  12 
but  N can  be  as  large  as  256  in  the  current  version  of  the  program , 

LAYR  is  the  number  of  the  layer  boundary  where  the  field  points  are 
located , And  it  is  the  only  one  of  the  variables  which  may  be  altered 
on  return  to  the  calling  program.  This  is  the  only  subroutine  in 
this  program  which  relies  on  memory  for  quantities  calculated  during 
previous  calls  of  it. 

Ths  logic  of  this  routine  is  not  complicated  when  it  is  broken  down 
into  the  parts  headsd  by  comment  cards.  The  first  part  evaluates 
the  line  Integral  of  B given  by  (40).  This  line  Integral  is  taken 
along  a straight  line  parallel  to  the  direction  of  incidence  to  the 
point  (p,  (j>)  in  question.  A given  path  of  integration  intercepts  a 
given  radius  p at  two  values  of  <|>  which  are  symmetric  about  90° , so 
the  range  of  covered  by  "do  6"  is  90°.  The  line  integrals  at  a 
given  value  of  <j>  are  calculated  in  implicit  loops  over  the  radial 
grid  points.  Between  iLatemente  2 and  3 the  line  integral  for  <j>  < 90° 
is  accumulated  in  the  variable  A as  p goes  through  all  radial  grid 
points  from  tw  down  to  the  one  in  question.  The  variable  B accumu- 
lates the  line  integral  from  this  point  down  to  but  not  Including  the 
point  of  closest  approach,  which  is  represented  by  the  variable  RS. 

(One  will  recognise  that  ths  equations  for  A and  B have  the  form  of 
the  line  integral  in  polar  coordinates  whan  the  radial  coordinate  is 
the  variable  of  Integration.)  The  contribution  up  to  the  point  of 
closest  approach  is  added  after  statement  number  5.  Statement  number 
6 takes  account  of  the  fact  that  the  Integral  for  <t>>  90°  la  that  for 
<t><  90°  plus  twice  the  integral  from  that  point  to  the  point  of  closest 
approach.  Multiplication  by  CAYR  (krw)  removes  the  normalisation  by 
rw,  and  the  fact  that  BOK  is  B/(k  sin  a)  Implies  that  the  array  BDS 
la  /B  da  for  the  grid  of  angles  from  0 to  180°, 

In  the  volume  Integral  in  (21)  we  need  the  effective  wavenumber  kB, 
which  is  the  gradient  of  the  phase  angle  of  E0*Gj.  This  derivative 
is  obtained  numerically  in  terms  of  the  values  of  this  scalar  product 
at  five  points  in  the  radial  and  angular  grids.  Thus  the  point  in 
question  in  the  integration  for  (21)  must  always  be  retardad  by  one  in 
both  directions  as  we  iterate  p and  or  else  we  would  not  have  the 
required  five  points  about  the  point  in  question.  Thus  LM  is  the 
layer  boundary  in  question  and  if  it  is  not  inside  of  MINI,  the  inner- 
most point  at  which  the  fields  were  calculated,  NM  is  set  equal  to 
one  greater  than  the  number  of  angles  needed  to  cover  the  region  0 to 
180°  in  the  average  over  angles,  In  the  "do  30"  loop  over  these 
angles  we  usually  (unless  it  is  logically  superfluous)  start  by 
evaluating  the  scalar  products  Btt • G0  for  five  polarisation  combinations. 
The  arrays  U and  V are  tho  real  and  imaginary  parts  of  these  products 
and  their  middle  subscript  denotes  the  polarisation  in  the  order 
vertical,  horizontal,  linear  orthogonal,  circular,  and  circular  ortho- 
gonal. The  last  subscript  pertains  to  the  angle  grid  and  the  first 
is  designed  to  Bave  the  value  for  the  current  radial  gtid  point,  the 
previous  one  and  the  one  before  that  in  locations  L,  LA,  end  LB, 
respectively.  To  interpret  further,  the  ten  values  of  3,  after  comple- 
tion of  the  "do  9"  loop,  are  the  ten  possible  dot  products  between 
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the  real  and  Imaginary  parts  of  the  fields  induced  by  each  of  the  two 
principal  plane  incident  wave  polarisations.  Refer  to  the  definition 
of  the  12  components  of  E give  in  h)  Subroutine  STORF,  which  applies 
also  to  the  components  of  the  W array  in  this  routine,  to  decode  these 

equations.  In  the  "do  25"  loop  UA  and  VA  are  the  real  and  imaginary 

parts  of  EU*G0  at  the  point  in  question  for  evaluation  of  the  integral 
over  <{>,  The  variable  T is  the  square  of  the  magnitude  of  this  quantity. 

The  effective  wavenumber  ke  is  the  gradient  of  the  argument  of 

E0-G0  - (U  + i V)  e2ikz  C0B  « (94) 

where  U and  V have  the  meaning  of  the  corresponding  FORTRAN  variables. 

Elementary  analysis  gives 


ke*  - (2k  cos  a)2  + [(U  - V |U)2  + (U 

- V |^2/p2]/lE8-G0|4  (95) 

DUR  and  DVR  represent  rw  timeB  the  partial  derivatives  with  respect  to 
pj  DUP  and  DVT  represent  rw/p  times  the  partial  derivatives  with  re- 
spect to  4>.  Note  that  the  radial  grid  is  nonuniform,  so  generally 
three-point  numerical  differentiation,  except  when  only  two  points  are 
had  at  the  edges,  is  used  for  the  radial  variation.  The  angular  grid 
being  uniform,  the  two  points  surrounding  the  point  in  question  are 
used  and  A<j>  “ 4n/N.  As  we  Indicated  above,  the  contribution  of  the 
critical  point  (LM  - KRIT)  is  ignored  if  |E0*G„|  is  possibly  resonant 
<T>1). 

When  "accumulate  averages"  is  reached,  A represents  j E0  *G0 1 2 times  the 
nonconstant  part  of  S(ke).  The  variable  FGS  accumulates  that  quantity's 
integral  over  (j)  by  the  trapezoidal  rule  and  the  variable  DEP  accumulates 
the  same  quantity  weighted  by  the  depolarization  by  second  Born.  Multi- 
plication by  F st  completion  of  the  integrals  normalizes  them  by  S(2k), 
the  free-space  spectrum  function.  Both  are  printed  out  layer  by  layer, 
so  one  can  see  the  relative  magnitude  of  the  second  Born  depolarization. 
Then  DEP  is  added  to  FGS  using  the  logic  that  the  depolarization  in  the 
orthogonally  polarized  wave  scattering  adds  to  the  principally  polarized 
wave  scattering  and  vice  versa. 

Finally,  if  the  layer  boundary  called  is  the  outermost  (LAYR  - LAYERS), 
we  add  one  more  layer  and  run  through  almost  the  entire  routine  again 
because  of  the  lag  in  the  layer  boundary  of  interest. 

1)  Subroutine  RINT(JM)  is  designed  to  calculate  the  normalized  moments 
of  the  scattering  per  unit  length  of  wake,  by  performing  the  integra- 
tion in  (21)  over  the  radial  coordinate.  The  argument  JM  is  the 
number  of  the  innermost  layer  boundary  at  which  the  Integral  over  ^ 
has  been  completed  by  subroutine  PRINT. 

The  output  of  this  routine  is  a set  of  20  integrals,  defined  as 
follows: 


. -w,  Wtikuk'  .UwVWittiUV  v* 
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(96) 


un  N1*  W 
1 + (u/w2) 


d(p/rw) 


where 

2 IT 

{41}  s 1 / | E0  *G  | 2 S(ke)  (1  + d)  d<f>  (97) 

2 tt  S ( 2 k ) 0 

Ne  is  the  axial  value  of  the  mean  electron  density  and,  in  the  defini- 
tion of  {<)>},  d is  symbolic  of  the  depolarization  by  second  Born.  The 
precise  definition  of  (97)  is  such  that  | E. »G0 | 2 multiplying  d is  for 
the  opposite  polarization  to  that  multiplying  unity  (See  paragraph  k) 
Subroutine  AVGS).  The  quantity  {$},  which  is  the  Input,  is  repre- 
sented by  the  array  FGS.  The  output  SO  takes  un  ■ 1 in  (96),  SI  takes 
v/v0,  S2  takes  v,z/v02,  and  S12  takes  (v/vB)2,  The  quantity  v is  the 
mean  velocity  profile,  v'  is  the  velocity  fluctuation,  and  vB  is  the 
axis  velocity.  Each  output  has  five  values,  one  for  each  of  the  five 
polarization  combinations  carried  through  from  AVGS.  Before  returning 
to  the  calling  routine  wo  print  out  these  20  integrals  for  each  of 
two  normalizations.  The  first  represents  normalization  by  the  first 
Born  approximation  for  principal  polarization,  in  terms  of  the  quanti- 
ties BN,  CN,  DN  and  EN;  i.e>,  these  normalizing  factors  are  the  results 
of  (96)  if  {4}  ■ 1.  Since  it  is  interesting  to  examine  refraction 
effects,  the  20  integrals  are  also  printed  out  in  terms  of  their  nor- 
malization by  the  values  of  the  first  Born  at  critical  density  for  re- 
frarcion.  (Recall  that  in  MAIN,  AN  is  [Ne/(NC  sin  2a))2.) 

m)  Subroutine  ZINT(NZ,V)  performs  the  convolutions  defined  by  (22),  where 
NZ  is  the  number  of  axial  stations  at  which  input  values  of  dSn/dz 
exist  and  V is  cos  a. 

The  variables  DSDZ,  DTDZ,  and  DVDZ  are  the  input  values  of  dSn/dz  for 
n ■ 0,  1,  and  2,  respectively.  Referring  back  to  MAIN, whence  these 
quantities  come,  one.  can  see  that  essentially  all  normalizations  have 
been  removed.  DSDZ  is  radar  cross  section  per  unit  length  of  wake  in 
units  of  meters.  The  other  two  quantities,  when  divided  by  DSDZ,  are 
the  mean  and  mean  square,  dopplcr  velocity  at  the  axial  station  in 
question,  in  units  of  the  body  velocity  and  the  square  of  the  body 
velocity. 

The  logic  of  this  routine  is  based  on  the  various  situations  which 
can  arise.  Thus  there  is  a special  form  of  result  for  normal  inci- 
dence (|  cos  a|<10“5),  as  follows! 

sn<°>  "-Z  ^ dz  (98) 

where  only  one  output  axial  station  is  required  since  the  backBcattered 
pulse  shape  is  identical  with  the  transmitted  pulse.  If  only  one 
axial  station  is  input,  including  the  case  of  normal  incidence,  then 
the  output  is  simply  set  equal  to  the  input  since  the  convolution  is 
meaningless.  In  all  other  situations  the  convolution  is  represented 
as  follows: 


(99) 
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where  a is  the  range  resolution  (distance  between  half-power  points 
for  the  pulse  backscattered  from  a point  target)  divided  by  cos  a. 

Note  that  the  pulse  is  assumed  to  have  the  shape  of  the  cosine 
squared  function.  This  integral  is  evaluated  at  a set  of  uniformly 
spaced  values  of  z regardless  of  whether  or  not  the  input  is  on  a 
uniform  grid.  The  number  of  points  used  to  calculate  it  by  the 
trapezoidal  rule  depends  on  the  relative  widths  of  the  scattering 
function  (the  distance  between  the  first  and  last  input  stations) 
and  the  pulse  (between  2ero  points).  If  the  pulse  width  dominates 
we  use  the  larger  of  20  or  the  maximum  allowable  number  of  axial 
stations.  If  the  scattering  function  dominates  we  use  an  odd  number 
which  approaches  unity  as  the  relative  width  of  the  pulae  approaches 
zero  and  which  Is  Identical  with  the  number  used  in  the  others  cnee 
when  the  pulse  width  approaches  the  width  of  the  scattering  function. 

The  limits  of  the  range  of  z for  the  output  generally  cover  the  points 
where  the  result  goes  to  zero,  where  it  ie  assumed  that  the  input 
scattering  function  drops  to  zero  at  either  end  of  lta  range  of 
etatione. 

Interpolation  of  the  scattering  function  dSn/dz  is  required  in  general 
because  the  grid  of  input  will  not  necessarily  match  the  integration 
grid,  The  form  of  the  Interpolation  has  been  designed  to  try  to 
account  for  situations  where  the  function  varies  either  linearly  or 
exponentially.  If  either  Interpolated  value  is  zero  or  if  they  are 
within  an  order  of  magnitude  of  each  other,  linear  interpolation  le 
used.  Otherwise  the  interpolation  is  exponential;  i.e.,  It  is  linear 
with  respect  to  the  logarithm. 

C.  Input  and  Output 

The  Hating  for  MAIN  gives  a complete  definition  of  the  input  requirements, 
bo  very  little  need  be  added.  A title  card,  which  can  contain  any  symbols  one 
wishes  in  all  columns,  must  precede  each  set  of  input  data.  One  will  notice 
that  certain  conventions  are  used  in  naming  the  input  variables.  All  dimensional 
quantities  contain  a reminder  of  their  units  at  the  end  of  their  name.  With  the 
exception  of  FMHZ,  all  names  which  contain  Z refer  to  axially  varying  quanti- 
ties and  ell  but  NZ  are  arrays.  MUTZ  is  defined  as  (6ip  - 9)  in  the  spectrum 
function;  thus  iJj  can  take  any  of  the  values  5/3,  11/6,  2,  13/6,  7/3,  and  5/2. 

The  program  will  return  to  the  beginning  of  MAIN  and  look  for  new  data  if,  after 
completion  of  all  calculations  for  a given  case,  the  number  of  cases  bo  far 
done  in  this  run  is  less  than  whatever  value  KASKS  has  been  given.  In  the 
Honeywell  6060  one  need  not  repeat  namelist  inputs  which  are  no  different  from 
what  is  currently  in  memory,  so  KASES  needs  to  be  input  only  with  the  first  set 
of  data.  In  the  definition  of  POLAD,  "plane  of  Incidence"  1b  the  plane  contain- 
ing both  the  wake  axle  and  the  radar  line  of  sight.  All  profile  functions, 
which  are  denoted  by  PRO  at  the  end  of  the  name,  are  defined  as  the  given  quan- 
tity at  some  radius  divided  by  a normalizing  quantity.  For  both  EMPRO  and  EPPRO 
the  normalizing  quantity  is  the  axial  value  of  the  mean  electron  density.  The 
collision  frequency  and  velocity  profiles  are  normalized  by  the  axial  value  of 
collision  frequency  and  velocity,  respectively.  In  general,  however,  it  does 
no  harm  to  use  inconsistent  normalizations;  e.g.,  one  could  input  a set  of 
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values  o£  IMPRO  which  represents  a profile  which  1b  different  from  unity  at 
the  axis  with  no  adverse  effect  on  program  operation.  The  EPPRO  often  will 
differ  from  unity  at  the  axis  if  N'  j*  Ne.  The  methods  of  inputting  various 
types  of  profile  functions  are  detailed  on  the  listing.  Note  that  when  a 
smooth  profile  function  is  to  be  Input,  the  data  correspond  to  evenly  spaced 
radial  points  starting  with  the  axis  and  ending  with  the  outer  edge.  When 
using  this  option  it  is  best  to  go  to  the  limit  and  input  eleven  points  in 
the  profile.  Obviously  the  meaning  of  the  variable  VVBZ  is  arbitrary  and  some 
other  velocity  than  that  of  the  re-entry  body  could  be  used  as  a normalizing 
factor. 

Other  than  those  already  mentioned , there  are  relatively  few  restrictions 
on  the  values  of  the  input  quantities.  ASPECD  must  not  be  zero  (or  such  that 
its  sine  is  zero).  CZPSEC,  EMZPCC,  PMHZ,  RESM,  and  RZM  must  be  greater  than 
zero.  The  values  of  EPPRO  may  not  be  such  that  the  profile  of  electron  density 
fluctuation  is  identically  zero.  Finally,  the  array  ZM  must  have  unequal,  In- 
creasing values. 

Appendix  B is  a copy  of  the  complete  output  for  a sample  CQBe,  which  can 
serve  as  a user’s  test  case.  The  output  starts  on  a new  page.  The  top  line 
gives  the  case  number  In  this  run  and  the  title  card  characters.  The  format 
for  this  line  is  given  in  statement  number  5000  of  MAIN.  We  have  used  a title 
card  spanning  all  80  columns,  to  illustrate  where  they  appear  on  this  line.  All 
data  for  the  input  variables  are  then  written  in  NAMELIST  format. 

The  next  new  page  gives  the  data  used  for  the  firBt  axial  wake  station. 

The  meanings  of  most  of  the  items  should  be  obviouB,  Each  "changed  to"  entry 
refers  to  the  item  to  its  left,  which  has  been  adjusted  by  PROFS  to  avoid  spur- 
ious effects  of  nearly  critical  density  in  a layer  (format,  including  column 
headings:  1000,  MAIN).  In  the  profile  functions  the  relative  radius  refers  to 
the  layer  boundaries,  which  are  called  "mesh  pt,"  In  the  listing.  The  next 
three  columns  fall  under  the  heading  of  electron  density.  All  other  columns 
are  self-explanatory  (format  for  each  row  of  data:  2000,  MAIN).  Note  that 
various  possible  kinds  of  input  profiles  are  illustrated  by  the  inputs  used  for 
this  case.  The  mean  electron  density  is  approximately  equal  to  a Gaussian  and 
the  rms  fluctuation  is  a step  function  to  half  the  wake  radius,  The.  collision 
frequency  is  uniform  and  the  two  velocity  profiles  linearly  decrease  by  a 
factor  of  ten. 

The  next  page  of  output  gives  the  summary  of  the  results  of  the  calcula- 
tion of  em  fields  in  the  mean  wake  (format  for  heading:  3000,  MAIN).  The  first 
column  shows  the  mode  orders  used.  The  next  eight  columns  give  each  of  the  two 
scattering  coefficients  for  each  incident  wave  polarization.  In  each  pair  of 
columns  the  first  is  the  real  part  and  the  second,  the  imaginary.  These  scat- 
tering coefficients  are  defined  in  the  discussion  on  subroutine  BCO  concerning 
equation  (89).  The  last  column  shows  that  the  wave  propagation  was  not  cut 
off  in  this  case  (format  for  each  row  of  data:  1000,  BCO). 

The  next  page  gives  the  far  field  coherent  scattering  pattern,  where 
"azimuthal  angle"  refers  to  <j)  and  the  firBt  row  of  the  table  is  for  the  forward 
scattering  (format  for  headings:  2000,  PHINI;  format  for  rows  of  numbers:  1000, 
FFSUM) , The  complex  forward  scattered  fields,  with  the  phase  relative  to  the 
incident  wave,  are  written  at  the  bottom  of  the  table  (format:  2000,  FFSUM). 
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The  next  page  is  the  output  of  subroutine  AVGS.  It  shows  the  average 
over  <(>  of  |E0  »G . 1 2S(ke)/S(2k)  for  each  polarisation  combination.  The  five 
columns  of  depolarization  factors  include  the  second  Born  in  this  average  as 
well  (format  for  headings:  4000,  PHINT;  format  for  data:  1000,  AVGS).  The 
data  for  this  sample  case  show  a very  interesting  phenomenon  on  this  page. 

The  critical  point  in  the  mean  electron  density  profile  is  at  point  # 33, 
where  AVGS  has  set  all  results  equal  to  zero  to  avoid  the  resonance  effect. 
However  this  wake  ijs  in  the  diffraction  regime,  which  permits  wave  propaga- 
tion without  cutoff  all  the  way  to  the  axis.  Diffraction  also  spreads  the 
resonance  effect  over  a band  of  radii  about  critical.  But,  as  can  be  seen 
from  the  data  either  Bide  of  critical,  the  phase  varies  ao  rapidly  in  this 
ragion  that  the  spectrum  function  drops  off  as  the  critical  point  is  approached. 

The  next  page  gives  the  lowest  three  moments  of  the  scattering  spectrum 
(format  for  page  headings:  1000,  RINT).  In  the  two  normalizations  by  the  Born 
approximation,  even  the  cross  polarized  results  art  divided  by  the  first  Born 
approximation  for  direct  polarization.  Also,  "first  squared"  refers  to  the 
contribution  to  the  second  moment  of  the  profile  of  the  square  of  the  mean 
velocity  (formats:  2000  and  3000,  RINT) , The  heeding  "absolute"  means,  for 
the  zeroth  moment,  radar  crosa  section  par  unit  length  of  wake  in  meters.  For 
the  first  or  second  moment  the  unite  are  matava  times  the  body  velocity  or  its 
square,  respectively  (format:  6000,  MAIN).  In  "dopplar/body  velocity",  the 
term  "mean"  is  aalf  explanatory,  "Spread"  is  defined  ae  the  standard  devia- 
tion of  the  dopplar  spectrum  in  units  of  the  body  velocity;  i.e., 

"spread"  5 (100) 

The  above  five  pages  of  output  are  repeated  for  each  ajrial  station  which 
waa  input.  Note  that  In  the  sample  run  wa  have  illustrated  a three-step  func- 
tion redial  profile  In  the  second  axial  station.  Also,  note  that  in  the  second 
and  third  stations  we  have  used  electron  denoitles  which  decrease  by  four 
orders  of  magnitude.  The  result  it'  show  that,  indeed,  the  Born  approximation  is 
approached  in  this  limit. 

The  final  pages  of  oucput  give  the  backBcattered  pulse  shape  functions 
for  each  of  the  five  polarizations  (format  for  headings:  1000,  ZINT).  The 
pulse  shapes  are  expressed  in  terms  of  radar  cross  section  (dBstn)  and  doppler 
mean  and  spread  (as  defined  above)  versus  apparent  axial  station  (format  2000, 
ZINT).  A page  is  turned  after  the  first  13  stations  and  then  after  each  14 
more  (format  3000,  ZINT), 

D.  Summary  of  Checkout  Procedures  Used 

The  reliability  of  this  computer  code  for  doing  the  calculations  for 
which  it  waa  designed  has  been  proven  by  a series  of  checkout  runo  during  Its 
development,  The  checkout  waa  done  at  various  levels  in  the  program,  starting 
with  subroutines  which  do  not  cull  other  subroutines  and  working  up  to  tha 
entire  program. 

The  Lagrange  interpolation  in  subroutine  FIT  was  checked  by  exercising  it 
on  two  profile  functions,  a Gaussian  and  a Gaussian  multiplied  by  a quadratic. 
The  input  data  were  taken  from  the  appropriate  pointB  on  the  true  underlying 
function.  The  number  of  points  used  was  varied  from  on*,  through  eleven  and 
the  accuracy  and  reasonableness  of  the  output  were  verified  by  Calcomp  plots 
as  well  as  numerical  printout.  For  each  of  these  functions,  one  can  not 
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visually  detect  any  difference  from  the  true  curve,  on  a graph  scale  of  about 
13  cm  height,  as  long  as  eight  or  more  points  are  input.  The  sample  output 
discussed  above  verifies  that  the  step-function  type  of  output  is  correct  also. 

* 

Subroutine  PROFS  was  checked  out  by  hand  calculation  of  its  outputs  for 
a few  selected  caseB. 

Subroutine  HANK  was  checked  out  by  comparing  its  output  with  tabulated 
values  for  a number  of  arguments  ranging  from  0.2  to  50.  The  results  were 
accurate  to  within  one  part  in  at  least  10”  or  within  an  absolute  value  of  10“® 
or  lens  if  the  result  is  much  less  than  10~6.  In  the  process  of  these  checkout 
runs  we  developed  the  convergence  criterion  which  is  built  into  this  routine. 

These  same  statements  apply  to  subroutine  BIN  as  well,  except  that  imaginary 
and  complex  arguments  also  were  used. 

Subroutine  PROP  waa  checked  out  by  comparing  itB  results  with  the  direct 
calculation  of  the  cross-products  of  Bessel  functions,  pn,  qn,  rn,  and  sn  in 
9,1.32  on  page  361  of  Abramowitz  and  Stegun  (op,  cit.).  (C  - - Trt0rn/2, 

D “ ^SoPn/2*  C'  - - TT£0sn/2,  D'  ■ irC0qn/2,  a - £0,  b ■ f„)  Subroutine  HANK 
was  the  standard  of  comparison  for  real  arguments;  library  routines  for  modi- 
fied Bessel  functions  were  the  standard  for  imaginary  arguments;  and  selected 
tabulated  values  were  the  standard  for  certain  real,  imaginary,  and  complex 
(using  Kelvin  functions)  arguments.  Initially  the  eerier*  expansions  of  the 
matrix  elements  in  PROF  were  formulated  in  terms  of  power  series  in  Y with  the 
coefficient  of  calculated  by  recurrence  relations.  This  approach  began  to 
lose  accuracy  when  |£0|  exceeded  10  or  20,  and  it  encountered  exponent  overflows 
in  the  Honeywell  6060  for  arguments  greater  than  50,  The  formulation  described 
above,  in  which  the  power  of  Y is  included  in  the  coefficient  which  is  calcu- 
lated by  recursion,  significantly  improved  the  performance.  Now,  as  Btated 
above,  we  can  go  to  values  of  [ | | (1  + Y) ] up  to  88,  which  is  the  ultimate 
single  precision  limit  of  the  Hotiewell  6060,  with  accuracy  generally  of  lO"'® 
or  better.  Values  of  Y UBed  in  these  runs  were  0.02,  0.05,  0.1,  0.2,  and  0.3. 

An  earlier  series  of  experimental  runs  using  a more  primitive  formulation  had 
shown  that  the  series  would  not  always  converge  if  Y waa  greater  than  0.3,  how- 
ever this  conclusion  may  now  be  conservative. 

Subroutine  WFP  and  the  routines  it  calls  (HANK,  BIN,  PROP,  BCD,  and  STORF)  were 
checked  out  by  doing  a sot  of  calculations  for  a uniform  cylinder,  'ilia  stan- 
dard of  comparison  waa  a specially  written  routine  using  Besnel  functions  in- 
side the  cylinder,  It  was  found  that,  even  for  zero  collision  frequency,  the 
accuracy  is  consistent  with  the  computer  word  size. 

Subroutine  FFBUM  was  checked  out  by  putting  in  the  Fourier  expansion  co- 
efficients for  the  incident  plane  wave  and  noting  that  itB  output  faithfully 
reproduced  the  plane  wave. 

Subroutine  PHINT,  which  calls  FFSUM  and  AVGS,  was  checked  by  doing  the  case 
of  zero  electron  density  In  subroutine  WPP  and  calling  PHINT.  This,  of  course, 
produces  the  incident  wave  as  the  output  of  FFSIJM.  AVGS  should  produce  renults 
which  are  unity  for  the  three  principal  polarizations  and  zero  for  the  ortho-  ; 

gonal,  This  result  was  obtained  except  that  zero  is  about  lO""^  for  linear 
orthogonal  and  10“®  for  circular  orthogonal. 

■( 

Subroutine  RINT  ves  checked  out  by  input  of  analytical  functions  for  which  | 

the  integrals  are  known  in  closed  form.  i 
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Subroutine  Z1NT  was  checked  out  by  putting  In  various  functions  of  z which 
can  be  convolved  with  the  pulse  shape  function  in  closed  form.  These  runs  in- 
eluded  cases  in  which  the  scattering  function  ranged  from  much  shorter  to  much 
longer  than  the  radar  pulse.  They  included  cases  in  which  the  scattering  func- 
tien  was  specified  in  terms  of  as  few  as  one  point  to  as  many  as  20.  Cases 
were  done  also  to  verify  the  program's  logic  for  normal  incidence. 

The  sample  run  discussed  above  serves  as  a checkout  for  the  entire  program. 

At  the  third  axial  station  the  electron  density  is  less  than  10"?  of  critical 
density  and  tha  principal  polarisation  rasults  are  1.0005  times  the  Bom  approxi- 
mation while  those  for  linear  and  circular  orthogonal  are  down  by  at  least  101® 
and  10^3,  respectively,  The  factor  1.0005  comae  from  the  numerical  differentia- 
tion to  obtain  the  affective  wavenumber,  which  therefore  does  not  exactly  equal  2k. 
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IV . SAMPLE  RESULTS 


A.  Comparisons  with  Experimental  Data 
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Shkarof aky  at  al.18  have  reported  experimental  data  on  microwave  scattering 
from  a turbulent  arc  jet.  Their  paper  Is  an  excellent  subject  for  analysis  by 
our  theory  since  it  gives  all  the  necessary  Information  on  the  plasma.  Radial 
and  axial  variations  of  mean  and  rms  electron  de-ity  are  given  in  terms  of  grapha 
of  the  profiles  at  nine  axial  stations,  and  the  turbulence  spectrum  function  also 
is  given.  Figure  4 ahowa  the  comparison  of  this  theory  with  Shkarof sky's  data. 
Mots  that  they  reported  the  absolute  backacatter  power  while  our  model  gives 
radar  cross  section  per  unit  length,  so  we  have  displaced  their  data  to  force  it 
to  agree  with  the  Born  approximation  at  low  alectron  density.  The  polarisation 
le  the  one  parameter  not  explicitly  stated  in  ref.  18;  but  even  in  the  worst 
case  that  the  polarization  was  parallel  to  the  plane  of  incidence,  the  agreement 
between  this  theory  and  the  data  is  excellent.  (Note  that  in  these  calculations 
fifty  angular  modes  were  used  on  a grid  of  256  angles.) 

Cross-polarized  backscatter  results  are  shown  in  Figura  5 in  terms  of  tha 
ratio  to  the  principal  polarization.  The  agreement  with  the  data  is  still  good, 
although  only  within  a factor  of  three.  The  moat  satisfying  aspect  of  both  sets 
of  results  is  the  very  good  agreement  in  the  shapes  of  the  curves  all  the  way 
through  critical  density.  Note  that  the  background  medium  contributes  to  depolar- 
ization also,  a fact  which  has  been  noted  by  Halseth  and  Sivaprasad. 11  The  dashed 
lines  in  Figure  5 show  that  this  contribution  is  quite  significant  in  this  case, 
being  about  equal  to  that  from  the  second  Born  approximation  when  the  plasma  is 
overdenae.  (Special  changes  were  made  in  the  program  to  output  these  results.) 

Graf  et  al.20  have  studied  the  effectB  of  incidence  angle  on  refraction  in 
a turbulant  laboratory  plasma.  Their  data  for  backscattering  at  40°  aspect  angle 
(between  the  line  of  incidence  and  the  plasma  axis)  follow  the  first  Born  square 
law  dependence  on  the  rms  electron  density.  These  data  points  extend  up  to  an 
rms  electron  density  of  0.7  times  critical  density.  They  do  not  show  backscatter 
data  for  20°  aspect  angle  although  they  state,  "At  an  aspect  angle  of  20°,  the 
backscatter  cross  sections  deviated  ...  from  the  slope-of-two  curve  ...  about 
3db  when  n'  was  0.15  nc."  Regrettably,  they  do  not  report  much  detail  about  the 
conditions  of  their  experiment,  especially  regarding  the  background  (average) 
electron  density  distribution,  which  is  so  important  in  this  theory.  We  have 
made  calculations  using  radially  Gaussian  electron  density  distributions,  where 
the  standard  deviation  of  the  GauaHian  is  2.1  cm  for  the  average  electron  den- 
sity and  2,5  cm  for  the  fluctuations.  The  results  of  calculations  using  this 
theory  are  shown  in  Table  2. 


18.  I.P.  Shkarofsky,  A.K.  Ghosh,  and  E.N.  Almcy,  "Direct  and  Cross-Po.larizpd 
Backscatter  Power  from  a Turbulent  Plaema",  Plasma  Physics,  Vol.  14,  pp. 
935-950,  October  1972. 

19.  M.W.  Halseth  and  K.  Sivaprasad,  "Depolarization  of  Singly  Scattered 
Radiation  in  a Turbulent  Plasma",  Phys,  of  Fluids,  Vol.  15,  No.  6,  pp. 
1164^-1166,  Juno  1972. 

20.  K.A.  Graf,  H.  Guthart,  and  D.G,  Douglao,  "Refraction  Effects  in  Turbulent 
Media",  Radio  Science,  Vol,  9,  Nos.  8,  9,  pp.  777-787,  August-September 
1974. 
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FIGURE  5. 


COMPARISON  OF  THIS  THEORY  WITH  SHKAROFSKY'S 
DATA  FOR  THE  RELATIVE  DEPOLARIZATION  OF  THE 
BACKSCATTERING . 
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The  parameter  N* /N  represents  the  axial  value  of  the  ratio  of  the  rtns  electron 
density  fluctuation  to  the  mean  electron  density,  and  It  can  be  seen  that  this 
Is  an  Important  parameter  in  the  theory.  Although  the  nominal  value  Is  given 
as  0.5,  the  theory  agrees  with  the  experiment  better  for  a larger  value. 

Graf  et  al.20  have  given  reiults  of  coherent  field  measurements  also. 

Figure  6 reproduces  their  rasul  for  meaauriSd  field  intensity  at  a distance  of 
156  cm  from  the  transmitter  when  the  plasma  axis  was  a perpendicular  bisector 
of  the  line  of  sight.  The  theory  gives  the  coherent  forward  scattered  field 
as  a byproduct,  which  falls  off  as  the  invern  square  root  of  the  far-field 
distance  because  of  tha  assumption  of  cylindrical  symmetry.  The  points  plotted 
on  the  graph  show  fchu  results, of  the  theory  for  the  total  field  intensity, 
assuming  that  the  Incident  field  falls  off  as  1/r.  The  square  symbol  denotes 
parallel  polarization  and  the  circle,  tranoverse.  The  open  symbols  show  the 
results  of  accounting  tor  only  the  mean  electron  density,  neglecting  the  effect 
of  scattering  attenuation.  The  filled  symbols  show  the  results  when  the  mean 
plasma  properties  nra  enhanced  to  account  for  scattering  attenuation,  assuming 
that  N'/N  « 0.5.  The  theory  gives  more  diffraction  than  the  data  show,  although 
the  trend  with  electron  density  is  similar.  It  would  appear  to  be  anomalous 
that  the  transmitted  field  is  increased  when  scattering  attenuation  is  accounted 
for}  but  this  result  is  caused  by  the  manner  in  which  the  calculation  is  done, 
in  that  quantities  are  added  to  the  mean  electron  density  and  collision  fre- 
quency. Apparently,  more  diffraction  than  attenuation  is  added. 

Attwood21,22  has  reported  the  results  of  backscattering  experiments  for  a 
wide  range  of  aspect  angles  and  electron  densities  either  side  of  critical. 

But  this  theory  comes  nowhere  near  explaining  the  extreme  degree  of  dependence 
of  his  data  on  these  parameters,  especially  at  and  above  critical  density.  This 
is  probably  caused  by  the  fact  that  his  plasma  does  not  fit  the  description 
aBBumed  heroin,  of  a random  medium  embedded  in  a steady,  cylindrical  background 
medium. 


21.  D.  Attwood,  "Microwave  Scattering  from  Underdense  and  Overdense  Turbulent 
Plasmas",  Physics  of  Fluids,  Vol.  15,  No.  5,  pp.  942-944,  May  1972. 

22.  D.  Attwood,  "Microwave  Scattering  from  an  Overdense  Turbulent  Plasma", 
Physics  of  Fluids,  Vol.  17,  No.  6,  pp.  1224-1231,  June  1974. 
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Axial  mean  electron  denaity/critlcal,  Ne/Nc 


FIGURE  6.  COMPARISON  OF  THIS  THEORY  WITH  GRAF’S  DATA  ON  THE 
COHERENT  FIELD  INTENSITY  TRANSMITTED  THROUGH  A 
TURBULENT  FLAME  AT  NORMAL  INCIDENCE. 


B.  Some  Parametric  Results 


In  the  previous  study  (see  Ref.  1)  the  importance  of  finite  wavelength 
effects  was  evaluated  In  some  detail.  These  effects  refer  to  many  of  the 
considerable  differences  which  are  found  when  the  results  of  rigorous  em  wave 
propagation  in  a cylindrical  plasma  are  compared  with  the  results  of  vay 
tracing.  Thus  that  study  identified  the  diffraction  effect,  which  was  so 
named  because  it  represents  the  departure  from  geometrical  optics  laws  as  the 
wavelength  is  increased.  Furthermore,  in  that  Btudy,  it  was  shown  that  the  j 

diffraction  regime  occurs  when  the  parameter  koe  sin  a is  small,  where  oe  is  ( 

the  width  (standard  deviation)  of  the  radial  profile  of  the  mean  electron  density 
distribution.  In  other  words  the  standard  of  comparison  for  the  wavelength  is 
the  width  of  the  electron  density  distribution  projected  onto  the  line  of  in- 
cidence. An  important  corollary  of  this  result  is  the  fact  that  low  radar 
aspect  angle  can  be  a sufficient  condition  for  diffraction,  Our  objective  in  ! 

this  subsection  is  to  determine  how  significantly  the  Improvements  embodied 
in  the  present  model  affect  the  results  both  in  the  diffraction  regime  and  in 
the  transition  between  diffraction  and  geometrical  optics. 

/ 

To  permit  direct  comparisons  we  have  used  the  same  radial  profile  functions  , 

as  for  most  of  the  prior  parametric  studies.  Thus  the  mean  electron  density 
is  a Gaussian  with  Pe  ■ 0.3rw;  i.e.,  exp  [-  p2/ (0.18rw2)) . The  rms  electron 
density  fluctuation  is  equal  to  the  mean  multiplied  by  [1  + p2/ (0.18rw2)] ; i.e.,  j 

it  has  an  Incipient  off-axis  peak.  The  mean  velocity  profile  is  Gaussian  with  ■ 

ov  ■ 0.94PB  a&,  The  rms  velocity  fluctuation  has  an  incipient  off-axis  peak  j 

and  the  standard  deviation  of  its  Gaussian  part  - 1,0787  ae.  We  use  a 
moderately  collisionless  condition  (\)/u  - 0.1)  and  a low  aspect  angle  (sin  a ■ i 

1/18).  The  scales  of  turbulence  are  kr^  - 0.021  and  kr0  ■ 0.21. 

Figures  7 through  10  show  the  RCS  (radar  crosa  section),  the  ratio  of 
orthogonal/principal  RCS,  the  mean  doppler  and  the  contribution  of  velocity 
fluctuations  to  tie  doppler  second  moment  as  functions  of  axial  electron  density  j 

for  circular  polarization , These  figures  are  all  for  kOe  sin  a « 0.03,  so  \ 

they  represent  the  diffraction  regime.  The  dashed  line  in  each  figure  repre- 
sents the  resultB  from  the  previous  study  and  we  conclude  that  the  improvements  j 

which  the  new  model  incorporates  (effective  wavenumber  calculation,  scattering  j 

attenuation,  and  depolarization  from  turbulence)  make  little  difference  in  the  I 

diffraction  regime.  1 

j 

Figures  11  through  14  show  the  same  four  quantities  as  functions  of  koe 
sin  a at  a value  of  Ne/(NC  sinza)  - 100.  The  only  significant  difference  be- 
tween the  old  and  the  new  results  is  the  polarization  ratio  for  kae  sin  a 
greater  than  about  2.  Thus  it  appears  that  in  the  diffraction  regime  the 
depolarization  by  the  meai.  wake  plasma  dominates  that  by  second  Born  scattering 
from  the  turbulence. 
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Normalized  axial  electron  density 


FIGURE  7 


RCS  VERSUS  ELECTRON  DENSITY  IN  THE  DIFFRACTION  REGIME 
THE  NORMALIZATION  IS  WITH  RESPECT  TO  THE  REFRACTION 
EFFECT;  I.E.,  THE  NORMALIZING  ELECTRON  DENSITY  IS 
N sin2a  AND  THE  NORMALIZING  RCS  IS  THE  BORN  APPROXIMA- 
TION AT  THIS  ELECTRON  DENSITY.  (FIGURES  7 THROUGH  14 
ARE  FOR  CIRCULAR  POLARIZATION  AND  THE  DASHED  CURVE  IS 
FROM  REFERENCE  1.) 
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Normalised  mean 
doppler  velocity  1 
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10”1  10°  101  102  103  10A 


Normalized  axial  electron  density, 
Ne/ (Ncein2a) 


FIGURE  9.  MEAN  DOPPLER  VELOCITY  IN  THE  DIFFRACTION  REGIME. 

THE  NORMALIZING  VELOCITY  IS  THE  BORN  APPROXIMfci 
TION;  X.B.,  THE  VELOCITY  PROFILE  WEIGHTED  BY  N71 
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Normalized  doppler 
second  moment  from  1 
velocity  fluctuations 

0 

10-1  10°  101  102  103 

Normalized  axial  electron  density, 
Ne/ (Ncein2«) 


FIGURE  10.  CONTRIBUTION  OF  VELOCITY  FLUCTUATIONS  TO  THE 
DOPPLER  SECOND  MOMENT  IN  THE  DIFFRACTION 
REGIME. 


kaeslnc* 


FIGURE  11. 


NORMALIZED  RCS  AS  A FUNCTION  OF  THE  ELECTRONIC 
WAKE  WIDTH  PROJECTED  ONTO  THE  LINE  OF  INCIDENCE 
AND  COMPARED  WITH  THE  WAVELENGTH.  (IN  FIGURES 
11  THROUGH  14  N /<NeBin2a)  - 100.) 
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Normalized 


FIGURE  13.  NORMALIZED  MEAN  DOPPLER  VELOCITY  VERSUS 
NORMALIZED  ELECTRONIC  WAKE  RADIUS. 
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FIGURE  14.  CONTRIBUTION  OF 
DOPPLER  SECOND  ! 
ELECTRONIC  WAKE 
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V.  CONCLUSION 


A sophisticated  computer  model  for  the  radar  backscattering  from  turbu- 
lent re-entry- induced  ionized  wakes  has  been  described  In  sufficient  detail  to 
permit  the  user  to  modify  the  program  as  desired,  This  model  uses  the  basic 
concept  that  propagation  in  the  background,  average  electron  density  distribu- 
tion determines  both  the  local  fields  and  the  Green’s  function  for  scattering; 
and  this  effect  is  formulated  in  a rigorous  and  general  way  by  assuming  that 
this  background  is  a cylindrical  medium.  At  this  stage  of  development,  however, 
the  application  of  this  model  to  the  turbulence  scattering  is  relatively  crude, 
being  of  the  type  generdly  used  in  so-called  distorted  wave  Born  theories. 

Thus  much. room  for  improvement  exists. 

Meny  Improvements  and  modifications  are  possible  without  basic  theory 
changes.  For  example,  one  could  generalize  the  program  for  blstatic  scattering. 
Thio  would  take  major  programming  changes  but  it  does  not  require  any  new 
theory.  Generalization  for  anisotropic  turbulence  would  be  relatively  straight- 
forward, and  allowance  for  various  types  of  radar  pulses  is  almost  trivial. 

The  same  can  be  said  for  Including  the  effects  of  antenna  patterns;  i.e.,  non- 
uniform  illumination. 

This  model  probably  could  also  be  applied  to  the  problem  of  radar  scatter- 
ing from  the  turbulent  plume  of  a rocket  exhaust.  The  main  constraint  would 
be  that  it  be  in  a regime  of  flow  that  simulates  the  turbulent  wake  flow.  At 
high  altitude,  where  the  plume  is  laminar  and  npherl.cal,  the  canonical  problem 
to  be  solved  would  be  scattering  from  spheroidal  media.  The  presence  of  strong, 
nonrandom  irregularities  in  the  plume,  such  as  those  created  by  shock  diamonds, 
would  require  yet  different  techniques  such  as,  perhaps,  the  geometric  theory 
of  diffraction. 

Various  basic  theory  improvements  are  possible  in  the  problem  of  turbulent, 
cylindrical  wake  scattering.  There  is  every  reason  to  expect  that,  with  the 
known  rigorous  solution  for  the  effect  of  the  background,  one  could  do  the 
statietlcs  of  the  scattering  from  the  random  component  of  the  medium  in  a rigor- 
ous way,  At  least  one  could  do  a Monte  Carlo  calculation.  Also,  the  quasl- 
optlcal  approximations  (and  their  accompanying  ad  hoc  assumptions)  borrowed 
from  the  distorted  wave  Born  approximation  could  ba  removed  by  carrying  out 
higher  order  iterations  in  the  Neumann  series.  Of  course  the  difficulties  of 
doing  this  are  not  trivial  because  one  would  have  to  get  near-field  solutions 
for  propagation  in  the  background  medium.  Obviously  this  is  easier  to  do  in 
the  diffraction  regime,  where  fewer  angular  modes  are  needed,  but  this  Is  where 
the  significance  of  this  model  is  the  greatest. 

In  conclusion,  this  model  provides  not  only  a valuable  tool  for  applications 
but  also  a framework  for  future  research  in  uhe  theory  of  scattering  from 
random  media. 
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